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RELAXATION  COMPUTATION  OF  TRANSONIC  FLOWS 
AROUND  WINGS  WITH  BLUNT  LEADINC-EDGE  AND  DISCUSSION 
ON  ITS  STABILITY  AND  CONVERGENCE 


Northwestern  Polytechnical  University  * 
j Zheng  Yuwen  and  Luo  Shijun 

ABSTRACT 

In  this  paper,  the  exact  velocity  potential  equation  and 
the  exact  boundary  conditions  were  used  around  the  blunt 
leading-edge  of  a  wing,  and  the  velocity  potential  equation  with 
small  perturbation  in  the  transverse  direction  and  large 
perturbation  in  the  longitudinal  direction  together 
with  its  corresponding  boundary  conditions  were  used  In  other 
areas  to  obtain  the  solution.  Numerical  example  1  was  for  a 
rectangular- wing  with  an  aspect  ratio  X=12,  airfoil  NACA0012, 
free  stream  Mach  number  M»=0.63,  and  attack  angle  a=2c.  The 
calculated  pressure  distribution  of  the  root  section  was  close 
to  the  exact  numerical  subsonic  solution  (Sells,  1968).  Example  2 
was  an  experimental  wing  NACA  RM  A51G31  having  airfoil  NACA  64A010 
which  is  perpendicular  to  the  1/4  chord  line  with  a  sweepback 
angle  x=3,  taper  ratio  q*2,  M«=0.4,  0.8,  0.9,  and  a=2°. 

The  computed  results  were  very  close  to  the  experimental  ones. 

In  this  paper,  we  establish  the  stability  conditions  of 
the  linear  relaxation  with  improving  iteration  of  the  transonic 
velocity  potential  difference  equation  with  small  steady  pertur¬ 
bation  under  the  assumption  of  local  linearization,  and  the 
conditions  for  the  convergence  of  the  relaxation  solution  to 
the  original  differential  equation  solution.  These  conditions 
more  or  less  agree  with  the  numerical  experiments. 

"Received  December,  1981. 
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INTRODUCTION 


The  blunt  leading-edge  is  a  singular  point  of  the  classical 
small  perturbation  equation.  There  are  numerous  treatment 
methods  in  the  literature.  In  Reference  [1,2],  the  leading-edge 
was  avoided,  i.e.,  the  leading  edge  was  not  taken  as  a  mesh 
point.  Thus,  the  arrangement  of  the  mesh  had  a  large  effect 
on  the  calculated  results.  Reference  [3],  in  the  treatment  of 
the  leading-edge,  used  the  following  equation: 

m±dx-yt( xT,  *) 

J  *L  dX 

where  x^  and  x,p  are  the  chord  direction  coordinates  of  the 
leading  edge  and  the  trailing  edge;  y±  are  the  vertical  coordinates 
of  the  top  and  bottom  surfaces  of  the  wing;  and  z  is  the  span 
direction  coordinate.  The  integral  on  the  left  is  calculated 
by  the  trapezoidal  equation  according  to  the  mesh  used.  Using 
the  above  equation,  it  is  possible  to  solve  for  — at  the 
leading  edge.  In  Reference  [4],  based  on  the  mesh  used, 

-±0.2,  was  obtained  at  the  leading  edge  using  an 
extrapolation  method.  This  type  of  treatment  is  equivalent 
to  the  sharpening  of  the  leading  edge  which  has  some  arbitrariness. 

In  the  first  part  of  this  paper,  the  leading  edge  of  the 
wing  was  taken  as  a  mesh  point.  Using  the  directional  derivative 
equation,  the  exact  boundary  condition  of  the  blunt  leading 
edge  was  Inserted  into  the  exact  velocity  potential  equation. 

In  other  areas,  the  small  transverse  perturbation  velocity  poten¬ 
tial  equation  and  the  small  transverse  pertubation  boundary 
conditions  were  used.  This  method  avoided  the  shortcomings 
of  the  above  described  method  which  is  also  easy  to  apply. 

The  results  of  several  numerical  tests  showed  that: 
whether  the  linear  relaxation  of  transonic  small  perturbation 
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potential  flow  converges  or  not  and  the  convergence  speed  (if  it 
converges) are  related  to  the  relaxation  factor 
Here  we  need  to  use  theoretical  analysis  guidance.  In  the 
second  part  of  this  paper,  under  the  assumption  of  local 
linearization,  the  stability  of  linear  relaxation  with  improving 
iteration  was  analysed  using  the  von  Neumann  method.  Furthermore, 
the  convergence  of  linear  relaxation  improving  iteration  was 
discussed  using  the  separation  of  variables  method  for  the 
corresponding  differential  equation. 

1.  RELATION  COMPUTATION  OF  TRANSONIC  POTENTIAL  FLOW  AROUND 
WINGS  WITH  BLUNT  LEADING  EDGE 

1.  Basic  Equations 

Let  us  choose  a  rectangular  coordinate  system  oxyz.  The 
x-axis  is  parallel  to  the  wing  chord  and  the  z-axis  is  parallel 
to  the  wing  span.  At  the  blunt  leading  edge,  the  exact  pertur¬ 
bation  vel  city  potential  equation  was  used: 

(a1— i4,)l,M  +  (a1— «*)«PW  +  (a‘  — u)*)5?,,  — 2uw«P„— 2iru'P„  — 0  (D 


where 


a>-al+-~2-±(ql~ui-Vl-wt) 

»  —  (J.cosa  +q>„  v  »■  g.sin  a  +  cpn  u>  *»<p, 

y  is  the  adiabatic  index;  q®  and  a®  are  the  velocity  and 
sonic  speed  of  the  free  stream;  and  (u,v,w)  and  a  are  the  local 
velocity  and  the  sonic  speed. 
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At  a  non-leading  edge  point,  assuming  that  velocity 
perturbation  components  in  the  y  and  z  directions  are  small 
and  that  in  the  x-direction  may  not  be  small,  then  equation  (1) 
can  be  simplified  as 


(  1  0 


where 


1  — Mi.ore’g— 


1  -M‘= 


l  -• 


Y  -  1 


<2- 


Micoea<P, - — j— Ml'Pi 


M  is  the  local  Mach  number. 


(2) 


(3) 


2.  Boundary  Conditions 

Let  us  choose  the  coordinate  plane  oxz  in  the  wing  plane. 
The  boundary  condition  of  the  wing  surface  with  the  exception 
of  the  leading  edge  can  be  simplified,  under  the  assumption 
of  small  transverse  perturbation,  as 


±  o,  i)~(g_cnsa  +<P,(x,  ±  o,  *  )>J ^--g.siiia 

Similarly,  the  conditions  on  the  free  vortex  behind  the  wing  can 
be  simplified  as 

*£*,  +  o,  z  -  o,  z)  (5) 

*(x,  +  o,  2)-v(x,-0,  *)«  <P(*,,  +  o,  *)-  <P(*r,  -  o,  /)  (6) 


The  exact  boundary  condition  at  the  blunt  leading  edge  is 

<?>.—  —  g.cos  a  cm  X 

where  x  -  the  sweepback  angle  of  the  leading  edge  of  the  wing, 
n  -  the  normal  direction  of  the  leading  edge  (Figure  1). 
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Inserting  equation  (7)  into  the  directional  derivative  equation 
of  the  leading  edge,  we  get 

<P,»T,sin  X  —  fl_co8  0cos*X 
<*>,«■  tp.cosX  +q. cos  a  sin  X  cos  X 

where  t  is  the  direction  of  the  leading-edge  (Figure  1). 


By  isolating  the  boundary 
condition  of  the  leading  edge 
of  the  root  of  the  wing,  equation 
(8)  can  be  transformed  into 

«p,*»  — g.cos  a 
<P,«  0 

Similarily,  the  boundary 
condition  of  the  leading  edge 
of  the  wing  tip  is  also 
equation  (9). 


The  boundary  condition  of  the  far  field  of  the  wing  can 
use  the  small  perturbation  condition  which  Is  shown  in  §1.6  of 
Reference  [6], 


3.  The  Difference  Equation 


Let  us  choose  the  sequential  symbols  for  the  mesh  nodal 
points  in  the  x,y,z  directions  to  be  i,j,k.  At  the  blunt 
leading-edge  of  the  wing,  the  first  order  partial  derivatives 
?x  and  9z  are  calculated  from  equation  (8)  in  which  9t  uses 
the  central  difference  scheme.  This  means  that  it  is  assumed 
that  at  the  leading  edge  it  is  subsonic: 

(10) 
(11) 

where  Ah-A*//sinx  ,  and  are  the  step  lengths. 


'  AhMAh(Ah.i  + Ah) 
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The  first  order  partial  derivative  q>  ,  similar  to<pt,  is 
calculated  using  the  central  difference  equation.  The  second 
order  partial  derivatives  can  be  calculated  as  follows: 


®  fw _ Vhh&ZlLliLi-) 

Ax,.,  V  A*,.,  ) 


Ax* 


,  T„^i.*Ay/.i-<1P.,/,*(Ay).i  +  Ay))-fT,.,.l.»Ay) 
”  A>MAy;(AyM  +  A>,) 


1 


Ax,., 


— (cP»l/.;.»—  T,|,.  ,,/.*) 


(12) 


where  <px  and  <pz  are  computed  based  on 
order  partial  derivatives,  similar  to 
the  central  difference  equation.  A  . 

V,  1 


equation  (8).  Other 
t ,  are  calculated  us 
and  Azk  are  the  step 


first 

ing 

lengths . 


With  the  exception  of  the  leading  edge,  the  small  transverse 
perturbation  velocity  potential  equation  (2)  is  used.  The 
Murman-Cole  mixed  difference  format  is  used  in  the  difference 
equation.  On  the  surface  of  the  wing  and  the  upper  surface 
j*jw+o  of  the  free  vortex: 


j*  (  _ __ 

\  A  y,m 


) 


(13) 


Similarily,  the  expression  for  the  lower  surface  can  be 
written. 
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4.  Linear  Relaxation  and  Iteration  Computation 


The  set  of  difference  equations  is  solved  using  the  linear 
relaxation  and  improved  iteration  method.  The  details  of  the 
method  will  be  shown  in  Section  II.  Numerical  example  1  is 
a  rectangular  wing  having  airfoil  NACA0012  with  an  aspect  ratio 
of  \=12.  Choose  a  mesh  which  is  31x13x33-  The  wing  chord  is 
equally  divided  into  20  blocks  and  the  half  wing  span  is  divided 
into  29  blocks.  It  takes  20  seconds  per  Iteration  on  the 
655  computer.  M®=0.63,  a=2°.  By  choosing  a  relaxation  factor 
(o=1.0,  where  the  number  of  iteration  reached  n=1045,  we  obtained 

m»xl*P Jr*.*- I «=  0 . 15 x  10*4  ,  which  is  noted  as  /A«p/ 

The  pressure  distribution  of  the  wing  root  section  obtained  is 
very  close  to  the  exact  numerical  subsonic  solution^-1-^  as  shown 
in  Figure  2. 


Figure  2.  Pressure  distribution  of  the  root  section 
of  a  rectangular  wing. 

Key:  1)  calculation  in  this  work. 


Example  2  is  the  experimental  wing  in  Reference  [8]. 

The  airfoil  perpendicular  to  the  1/4  chord  line  is  NACA  64A010 
Xi/1,*450>  x=3,  n= 2.  Choose  a  mesh  of  31x13x19.  There  are  15 
sections  along  the  half  wing  span.  Each  section  has  11  points 
It  takes  10  seconds  for  each  iteration  on  the  655  machine. 

The  computational  results  are  shown  in  Table  1.  The 
computed  results  agree  with  the  experimental  ones  as 
shown  in  Figures  3-5.  The  calculated  results  are  better  than 
those  reported  in  Reference  [4].  In  Reference  [4],  the  small 
perturbation  velocity  potential  equation  and  ■  -±0.2 

were  used  at  the  blunt  leading-edge. 


Table  1.  Computed  cases  of  NACA  RM  A51G31  Wing. 


a 

1 

« <M>  1 ) 

n 

IA*I 

0  4 

2* 

i 

1 

258 

0.7  *  10-* 

0.8 

2* 

0.9 

0.7 

481 

0.13X  10'* 

0  9 

2" 

0  9 

0.7 

667 

0.15X10-* 

8 


t-  v 


LO 


Figure  5.  Pressure  Distribution  of  NACA  RM  A51G31  wing 
M— 0.9,  o*2° . 

Key:  1)  computed;  2)  experimental  [7]. 


2.  DISCUSSION  ON  THE  STABILITY  AND  CONVERGENCE  OF  LINEAR  /8 

RELAXATION  AND  I L. PROVING  ITERATION 

1.  Discussion  on  Stability' 

In  order  to  analyze  the  stability  of  the  difference  equation 
of  the  linear  relaxation  improving  iteration  of  the  small  trans¬ 
verse  perturbation  velocity  potential  equation  (2),  it  is 

2 

assumed  that  the  coefficient  1*M  of  the  <p  xx  term  in  equation 
(2)  is  a  constant,  i.e.,  local  linearization.  Thus,  at  the 
local  subsonic  point  and  the  local  supersonic  point,  equation 
(2)  can  be  transformed  into 

o 
o 

By  choosing  the  relaxation  line  which  is  parallel  to  the  y-axis , 
the  difference  equations  of  linear  relaxation  and  iteration 
of  equations  (14)  and  (15)  are 


(14) 

(15) 
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Introduce  a  new  variable  t  which  corresponds  to  the  separation 
variable  n.  Then: 


(B) 


Let  the  step  length  of  t  be  At. 


leading  difference  form: 


a<P 


should  choose  the 


The  difference  equations  (19)  and  (20)  are  equivalent  to 
the  following  differential  equations: 

+  - 1  )  A«5  +(  ')u' “ 

The  above  equations  are  different  from  equations  (14)  and 
(15).  They  are  time  dependent  equations. 


(21) 


(22) 


(23) 


i  9 


How,  let  us  analyze  the  stability  of  equations  (19)  and 
(20).  Let  the  exact  solution  of  the  difference  equation  be  <p, 
and  the  numerical  solution  of  the  difference  equation  be  v+  5 
where  <5  is  the  error  introduced  by  the  numerical  computation. 

and  +  6  satisfy  the  same  difference  equation.  Because  the 
difference  equation  is  linear,  therefore,  6  also  satisfies  the 
same  difference  equation.  Based  on  the  superposition  principle 
of  linearity,  <S  can  be  decomposed  into  the  basic  solution 
with  an  expression  such  as 


b(x,y,  i,  i  (24) 

where  T  ,  0^,  02>  and  ^3  are  arbitrary  real  numbers  and 

<*0  is  a  function  of  these  numbers. 
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Substituting  equation  (24)  into  (19)  and  (20),  we  get 

(  2  1  Y  1  I  2  -tin*  I  1  ^  1  -COS0.A* 

!XaF  +  A/ m  2  +  A?/ - A? 

2  oj„t  0tAy  1  — cosp,A*  ,  ..  /sinP,Ax  sinfcAa  \ 
~a7~81"  2  Ax1  +  Ax’  ) 


1  +  2  sin»J^ 

.  A*1  Ay*  2 


.  1_  \  ,  1  -COs3,A* 

■  +  'SrJ+ - A?- 


+ A *.-^+ j-=^4£. +j>  ( jsei4£_+. 

Ay*  2  A*1  V  A*1 


sing.  A* 
A*1 


where 


"  -(^r- '  Xa?  +  a?  *  i0+ ^  «•  ^ 

S’(«“  IXijr  +  4?s,n'  ^  +  •Sr)-  4^r<x»8,a«i»'-^  +  ay  Si.' ■&&*. 


According  to  the  stability  condition  |e*,4,K  1 
(25)  we  obtain 


,  from  equation 


0  <a><  2 

For  equation  (26),  it  is  not  possible  to  find  an  u  to  satisfy 
the  stability  condition.  The  proof  is  shown  in  the  Appendix. 

2.  Discussion  on  Convergence 

Convergence  means  the  converging  of  the  solution  of  the 
time  dependent  differential  eauetions  (23)  and  (24)  to  the 
solution  of  the  steady  state  differential  equations  (14)  and 
(15).  Perform  variable  transformations  for  equations  (22)  and 
(23), 


,+  2  iT’+~T 


A* 
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We  obtain 


■&  )>•-(  +  <«> 
1  ( »••-[(  i  ■  * )i“ +  (~_ 'Wl*'- 0  (29) 

Equation  (28)  is  hyperbolic  and  equation  (29)  is  super  hyperbolic. 

By  using  the  variable  separation  method,  let 

<P  =F(  *  )C(  X,  y,  *)  (A) 

From  equation  (28),  we  get 


+((  g-J  1 ■  At&d&L 


(B) 


The  solution  is 


<F(*.  y,  *,  *)-=G„(x,  .y,  *)  +  2  O 

m  »  1 


where 


/>-* 


Af 

2_ 

A* 


Gm(x,y,z)  is  the  characteristic  function  of  the  boundary 
value  problem  of  the  following  equation 

G., -f  G^+G,, +  feLG  0 

.2 

^m  is  the  corresponding  characteristic  value.  Let 


(30) 


(3D 


(C) 


«<«<•"<&*  .<«•<*" 


CD) 


/10 
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^0  (x,y»z)  i3  the  solution  of  the  Laplace  equation 

Gb+Gjj+C,,*  0 

From  equation  (30),  we  know  that  when  0<w<2.  Hence,  the  convergence 
condition  of  a  local  subsonic  velocity  point  is  0<w<2.  The  converging 
speed  is  determined  by  the  minimum  value  of  the  real  part  of  the 
exponents  in  equation  (30),  i.e.  Re.  (p^).  From  equation  (31) 


where  Ax  and  Az+0 


lim^VCx,  y,  z,  r  )«<?,(*,  y,  g  ) 


We  can  see  that  along  with  increasing  uj,  p  increases  monotonically . 
Under  usual  conditions,  it  is  possible  to  choose  w  so  that  the 
square  root  in  the  expression  of  p  is  zero  in  order  to  raise  the 
converging  speed.  This  optimal  relaxation  factor  is 

- - - 

1 +7  t~t  (32) 

Similarly,  the  solution  to  equation  (29)  is  found  to  obtain 
the  convergence  condition  of  local  supersonic  velocity  points 


0  <®<  1  +' 


and  the  optimal  relaxation  factor 


1  n  *9 

Furthermore,  with  increasing  w,  the  converging  speed  increases 
monotoni cally .  When  Ax  =  Az-*-0,  the  optimal  relaxtion  factor  is 


I 


3.  Comparison  With  Numerical  Tests 


The  analysis  of  the  stabilization  of  the  above  local 
linearization  pointed  out  that:  For  three  dimensional  super¬ 
sonic  stream,  linear  relaxation  is  always  unstable.  However, 
many  numerical  computations  are  stable.  In  addition,  the  results 
analyzed  in  the  above  two  sections  agree  with  the  experience 
of  numerical  calculation. 

In  Reference  [3],  the  isolated  wing  of  experimental 
model  NASA  TND-83O  was  computed.  The  wing  was  triangular, 
x=60°,  the  flow  section  is  NACA  65  A  003,  M »  =1.05,  and 
0=2.2°.  For  various  relaxation  factors,  the  convergence 
conditions  are  shown  in  Table  2. 

lAVKlO** 

Table  2.  Convergence  cases  of  NASA  TN  D-83O  wing,  M«*1.05,  a=2 . 2 


KMAft* 

1 

1.0 

l.T 

I  1 

1.8 

l.T 

l.T 

1.0 

***««» 

2 

o.t 

#.7 

i  *» 

1.0 

1.0 

1.3 

1.5 

14* 

|  553 

|  231 

191 

194 

149  1 

70 

4 

|  465 

!  «» 

365 

339 

335 

219 

153 

Key:  Du  subsonic  velocity  point;  2)wof  supersonic  velocity 

point; 3)of  n  for  stopping  oscillation;  4)  n  for  reaching  the 
convergence  standard. 

When  u*1.6  for  a  supersonic  velocity  point,  the  computation 
is  divergent.  When  it  is  1.5,  convergence  is  the  fastest. 
This  agrees  with  the  convergence  analysis  in  Section  2. 


3.  CONCLUSIONS 

In  this  paper,  the  exact  velocity  potential  equation  was 
used  at  the  blunt  leading-edge  of  the  wing  to  overcome  the 


uncertainty  at  the  leading  edge  using  the  small  perturbation 
method.  The  pressure  distribution  of  the  wing  plane  obtained 
was  improved.  In  this  paper,  in  other  parts  of  the  flow  field, 
the  small  transverse  perturbation  velocity  potential  operation  was 
used.  It  saved  considerable  computer  time  as  compared  to  the 
use  of  the  exact  potential  flow  of  the  entire  flow  field.  It 
is  also  easier  to  be  applied  to  the  complex  wing-body  structure. 

In  this  paper,  under  the  assumption  of  local  linearization, 
the  stability  and  convergence  conditions  of  linear  relaxation 
and  iteration  were  established.  In  addition  to  individual 
conclusions,  it  agreed  with  the  numerical  tests.  It  is 
meaningful  in  the  sense  of  providing  some  guidance. 


APPENDIX  The  proof  of  the  constant  instability  of  linear  relaxa^ 
tion  iteration  in  the  local  supersonic  region. 


From  equation  (26) 

IBI'-UI1-  l(  “  ~  '  X  l-1  +  Si"'  S’!  V  +  A2i 

-  4<,  Sft  ♦  &  - *■&“)♦( 

j.  2 _ cosg,A^  V  1 

A*1  )  A*4 


A*1 
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(1)  When  u)*l,  choose  PtAy-PjA*-*  0  and 

When  O<0,A*<-|  W,  jBl1— 0 
Hence,  when  w*l,  it  is  unstable. 
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(2)  When  w<l,  choose  0»Ay-0,A*—  0  and  Ax—  Ax 

+  -^3in‘  J^asin’-^-cosP.Ax) 

when  O<0,Ax<-y~,  ]B|‘— 1^4|*<0 

Therefore,  when  w<l,  it  is  unstable. 

(3)  when  u»l,  let  e>o 

Choose  Pi  Ax Ay-  0  and  Ax--j^~ 

|Br-U!*-(  1  ~4  e  )~i  (  1  -COS0.AX) 

when  e>l/4,  \B\* -|yi|’<  0 . 

Choose  0,Ay— 0,Ax- 0  and  ax-a* 

|fl|»-|^|*-^r3in1  +(3  e  -  l  )cosp,Ax] 

When  e<l/4  and  B-^X  is  sufficiently  close  to  zero,  |B|*-|^|l<0 
Therefore,  when  e>0,  i.e.,  w>l,  it  is  unstable. 
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RELAXATION  COMPUTATION  OF  TRANSONIC  FLOWS 
AROUND  WINGS  WITH  BLUNT  LEADING-EDGE  AND 
DISCUSSION  ON  ITS  STABILITY  AND 
CONVERGENCE 

Zheng  Ywven  and  Luo  Shiftm 

( Northwestern  Polyiechnical  University) 

Abstract 

In  this  paper*  the  blunt  leading-edge  of  a  wing  is  taken  as  mesh  points, 
and  there  the  exact  velocity  potential  equation  with  central  difference  scheme 
and  the  exact  boundary  condition  are  used,  while  in  the  other  places,  the  ap¬ 
proximate  velocity  potential  equation,  which  assumes  small  perturbation  in 
the  transverse  plane  but  allows  large  perturbation  in  the  longitudinal  direction, 
and  the  corresponding  boundary  condition  are  employed. 

Two  numerical  examples  are  following! 

(  1  )  A  rectangular  wing  having  airfoil  NACA0012,  aspect  ratio  X  “12,  an¬ 
gle  of  attack  a“2',  free  stream  Mach  number  M»“0.63.  The  computed  pre¬ 
ssure  distribution  of  the  root  section  agrees  with  the  exact  numerical  subsonic 
solution  given  by  Sells  (1968). 

(  2  )  The  sweepback  wing  tested  by  NACA  RM  A51G31  having  airfoil  NA- 
CA64A010  which  is  perpendicular  to  1/4  chord  line  with  sweepback  angle 
X„4“45*,  X  “  3  and  taper  ratio  *1  “  2 ,  <*“2',  M„“0.4,  0.8  and  0.9.  The 
computed  pressure  distributions  agree  well  with  those  obtained  by  tests. 

Under  the  assumption  of  local  linearization,  the  stability  of  the  difference 
equation  in  line  relaxation  with  Seidel  iteration  is  studied  by  the  von  Neumann 
method  and  the  convergence  of  the  solution  of  the  differential  equation  equi¬ 
valent  to  the  above  difference  equation  to  the  solution  of  the  original  differen¬ 
tial  equation  is  discussed  by  the  method  of  separation  of  variables. 

The  following  conclusions  are  obtained! 

(  1  )  The  stability  condition  for  the  line  relaxation  with  Seidel  iteration  is 
0<w<2  at  locally  subsonic  points,  where  ®  is  the  relaxation  factor. 

(  2  )  At  locally  supersonic  points,  the  relaxation  is  always  unstable.  The 
convergence  conditions  are  as  follows.  Let  the  steps  Lx  (chordwise)  and  Lx 
(spanwise)  perpendicular  to  the  relaxation  line. 

(3)  0  <e><  2  ,  at  locally  subsonie  points. 

(O  o  <»<i  + — Ar 

,+(£- 

The  numerical  experiences  agree  with  the  conclusions  (  1  ),  (  3  )  and  (  4  ), 
but  do  not  agree  with  the  conclusion  (  2  ). 


rr~»  at  locally  supersonic  points. 
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DETERMINATION  OF  AERODYNAMIC  COEFFICIENTS  FOR  A 
RE-ENTRY  BODY  BY  MEANS  OF  AN  EXTENDED  KALMAN  FILTER 
AERODYNAMIC  RESEARCH  AND  DEVELOPMENT  CENTER  OF  CHINA 


Jiang  Quanwei,  Xu  Jinzhi,  and  Zhou  Shuying * 


ABSTRACT 

In  this  paper,  an  extended  Kalman  filter  method  was  used 
to  determine  the  major  aerodynamic  coefficients  of  re-entry 
bodies.  The  emphasis  was  placed  on  estimating  states  and 
parameters  using  the  measured  data  during  the  re-entry  flight 
under  the  condition  that  trajectory  observation  data  was 
absent . 

The  data  included  body  axial  angular  rate  and  acceleration 
obtained  from  the  rate  gyros  and  accelerameters .  A  mathematical 
model  was  established  based  on  six-degree-of- freedom  motion 
equations.  Both  ballistic  and  maneuvering  re-entries  were 
considered.  Numerical  simulation  and  actual  measurement 
conversion  showed  that  the  present  method  provided  more 
satisfactory  results. 


SYMBOLS 


Cm  resistance  coefficient 

Cm,  C*  derivatives  of  normal  force  and  lateral  force 

Cm„  C.i  derivatives  of  pitch  moment  and  yaw  moment 

Cmti,  Cm,  control  moment  derivatives 

Cm*  C„  damping  moment  derivatives 

C,m,  Cm,  roll  and  roll  damping  derivatives 

C,*,  C i>  bottom  diameter 

d  gravitational  acceleration 

0  flight  altitude 


*  Received  in  March,  1931. 
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rotational  moment  of  inertia 
atmospheric  model  parameter 
mass  of  the  re-entry  body 
accelerations 
angular  speeds 
reference  area 
flight  time 
velocity  components 
combined  speed 

distances  from  the  accelerometers  to  the  center  of  gravity 
pose  angles 


p  density  of  the  atmosphere 

p,  reference  density 

bQ,  t>r  inclination  angles  of  the  control  plane 
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1.  INTRODUCTION 

In  recent  years,  the  determination  of  aerodynamic  characteristics 
of  a  space  craft  from  measured  data  obtained  in  flight  using  varicv 
parametric  identification  methods  is  one  of  the  important  subject*, 
in  the  astronautical  industry.  This  work  has  important  signi¬ 
ficance  in  the  design  and  final  planning  processes  for  a 
spacecraft.  First  of  all,  it  is  based  on  the  measured  data 
obtained  in  an  actual  flight  environment.  Consequently,  the 
shortcomings  due  to  ground  equipment  such  as  insufficient  wind 
tunnel  simulation  can  be  remedied.  Secondly,  it  provides  real 
data  for  the  design  of  the  guidance  and  control  systems.  This 
is  especially  true  of  re-entry  spacecraft.  There  are  very  few 
experiments  and  the  cost  of  experiment  is  high.  Therefore, 
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it  is  extremely  important  to  obtain  the  most  possible 
useful  information  from  the  analysis  of  limited  data. 

Since  the  emergence  of  the  Kalman  filter  method^ 
it  has  been  widely  applied  in  aspects  such  as  communications 
and  control.  Especially  in  astronautical  engineering,  its 

r  2 1 

application  is  even  more  popular1-  .  However,  modern  Kalman 
filter  methods  were  only  used  to  determine  the  aerodynamic 
coefficients  of  a  flight  vehicle  in  recent  years.  In  Reference 
[3] ,  based  on  the  re-entry  vehicle  point  mass  differential 
equation  of  motion,  several  Kalman  filter  plans  were  presented. 
References  [4,5,6]  discussed  the  use  of  extended  Kalman 
filter  methods  to  determine  the  aerodynamic  characteristics  of 
tactical  aircraft.  Reference  [7]  presented  the  feasibility 
of  real  time  estimation  of  aerodynamic  coefficients  using  a 
Kalman  filter  method.  These  type  of  efforts  in  foreign 
countries  have  already  obtained  some  progress.  However,  they 
all  Included  trajectory  measurement  information  such  as  altitude, 
velocity,  and  position.  Their  work  did  not  involve  the  identifi¬ 
cation  of  the  aerodynamic  coefficient  of  a  re-entry  vehicle 
under  the  condition  that  trajectory  observation  data  was  absent. 

The  purpose  of  this  work  was  to  attempt  to  use  the  Kalman 
filter  method  to  solve  the  problem  of  determining  the  aerodynamic 
coefficients  of  an  re-entry  vehicle  in  the  absence  of  trajectory 
observation  data.  This  problem  has  a  background.  Experience 
showed  that,  due  to  various  reasons,  the  trajectory  observation 
data  of  re-entry  flight  tests  frequently  could  not  be  or  could 
only  be  partially  obtained.  If  an  estimation  method  could  be 
found  to  determine  the  important  aerodynamic  parameters  based 
on  the  acceleration  and  angular  rate  data  measured  on-board, 
it  would  be  alot  more  meaningful.  This  paper  is  the  theoretical 
simulation  and  the  actual  measurement  conversion  specifically 
with  respect  to  this  problem. 
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The  major  differences  between  this  paper  and  Reference  [7] 
are  as  follows: 

(1)  In  Reference  [7],  the  flight  altitude  was  already  known. 
The  H  in  this  paper  was  unknown  which  was  a  state  quantity. 

(2)  In  addition  to  trajectory  re-entry,  maneuvering  re-entry 
was  also  considered. 

(3)  The  method  in  this  paper  has  already  been  applied  in 
practice  while  Reference  [7]  was  a  feasibility  study. 


2.  EXTENDED  KALMAN  FILTERING  UNDER  EXPANSION  CONDITION 

The  Extended  Kalman  Filtering  method,  abbreviated  as  EKF,  ex¬ 
tends  the  use  of  linear  Kalman  filtering  to  the  non-linear  secondary 
filter  system.  With  regard  to  the  study  of  re-entry  flight,  the 
dynamic  equation  can  be  written  as 

y(  *  )-?c?(  t ),  £,  / )  (i) 

where  represents  the  state  vector  of  the  vehicle,  C  represents 
the  performance  parameters  such  as  the  aerodynamic  coefficients 
of  the  flight  vehicle,  and  £  is  the  non-linear  differentiable 
function.  During  the  period  of  consideration,  let  us  assume 
that  C  does  not  vary 


c -o 

Combine  equations  (1)  and  (2),  then  we  get 


IX  1  >' 

•?c?(  * ).'  r,  tr 

1 1 J 

2 

Define  a  new  expansion  state  vector  x  (t) 


(2) 


(3) 
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Also  assume  that 


/(*(  0.0 


'?££(<  ).  £,  O' 
o 


Then,  we  have  the  usual  state  equation 


*(0-/C*(0,  0  +  G(0?(0 


(O 


(5) 


where  W(t)  is  the  Gaussian  white  noise  whose  average  value 
is  0  and  spectrum  density  is  Q(t).  It  is  used  to  simulate 
process  noise. 

Let  us  assume  that  the  observation  vector  is  a  set  of  dis¬ 
crete  values  which  is  a  non-linear  differentiable  function  of 
the  state  vector 


+  !jc»  A-”*  1 »  2 , 

where  VR  is  the  measurement  noise,  which  is  expressed  by  the 
positive  state  random  vector  with  an  average  equal  to  0.  Its 
joint  square  difference  matrix  is  RR.  Then,  there  are  the 
following  extended  Kalman  filtering  iteration  equations: 


x(  «  )“»/(x(  O.  O 

p(O-FCxO),  OP(i)+P(<  )FrCx(  t  )>  0  +  C(  «)<?(  0<^(  O 
Kk- Pk{  -  )Hl  {HkPk(  -  m +Rk)-' 

Xk(  +  )**Xic(~)  +  AV{Zk— ?!icCXk(  — 

P«(  + )  -  <  /  -KkH* Cx«c - ))> p«(  - ) 


where 


FCxCO.  0> 


ttOLLh  o 

*£(<> 


X(»)«X(f) 


H«Cxx  (-)3 


X(/*)»X*(-) 


and  the  initial  condition  of  state  xQ  and  PQ. 


(7) 


(B) 

(C) 
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In  the  simulation  computation,  equations  (5)  and  (6)  are 
used  to  produce  the  observation  data. 


III.  STATE  MODEL  AND  OBSERVATION  MODEL 

The  mathematical  model  of  this  work  has  the  following 
characteristics: 

(1)  It  is  based  on  the  usual  body  axis  six-degree-of- 
freedom  equations  of  motion; 

(2)  The  measured  data  only  includes  the  acceleration  and 
angular  rate  data  on-board; 

(3)  The  re-entry  body  can  perform  maneuvering  flight, 
i.e.,the  aerodynamic  rudder  can  carryout  control; 

(4)  It  has  asymmetric  aerodynamic  characteristics; 

(5)  It  considers  the  effect  that  the  accelerometer 
is  not  located  at  the  center  of  gravity. 

The  mathematical  model  in  this  paper  has  the  following 
limitations: 

(1)  In  the  small  attack  angle  flight  trajectory  section, 
the  aerodynamic  coefficients  can  be  expressed  as  a  linear  function 
of  the  attack  angle; 

(2)  Atmospheric  density  varies  exponentially  with  altitude 
and  its  mathematical  model  is  already  known; 

(3)  There  is  no  systematic  error  in  the  measurement. 
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Corresponding  to  equation  (1),  the  actual  form  is 
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OS  N 

u  —  vr  —  mo  + - C 9  sin  4 

m 

i>**wp—ur+  Qi(cyt-y-+C,h,br  )  +  5 cos  4  sin  Y 
v»»ug-up--^-(c..-^-+C.^)+  0  cos  4  cosY 

Kx)MxMx>-+c-“69M1-x)i'  |  (8) 

'-(sr)W-rH»;)c‘'+c-!'  H 1  -Jr)M 

4*  4 cosY  —  r sin  Y 
4»(4sinY  +  rcosY  )/cos  4 
Y-  P  +tg  4  (  4 sin  Y  +  r  cosY  ) 

H»Ksin4  —  K-^-cos4oobY  -  »coa4sin  Y  ; 


The  auxiliary  relations  are 

P  =»P,exp(  —  W/fct) 

+  + 

0=0. 5P^* 

The  actual  form  corresponding  to  equation  (6)  is 

£C*(  t  )1  -  (»,»y»,p«r)x 

where 

«.  = +0^)— (/>*  + «*) +-^-«r 
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where  x^ ,  y2 ,  z2 ,  y^,  are  known  constants.  At  this  place, 
an  acceleration  sensor  is  installed  in  the  following  manner: 

the  axial  sensor  is  placed  on  the  longitudinal  axis  and  the  trans¬ 
verse  sensor  is  in  the  cross-section  passing  through  the  center 
of  gravity. 


IV.  RESULTS  OP  NUMERICAL  SIMULATION  AND  PRELIMINARY  REAL  DATA 
EXTRACTION 

The  numerical  simulation  corresponds  to  the  maneuvering 
re-entry  situation.  The  control  plane  moves  regularly  according 
to  6g=br=0.lsin(20 1  )  .  The  flight,  aerodynamic  and  observation 
characteristics  of  this  computational  example  are  shown  in 
Tables  l-1). 


Table  1.  Physical  Constants 
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9.81 
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Key:  1)  g  (m/sec);  2)pn  (Kgsec2^14 ) ;  3)  kn(m);  4)  X,  (m); 

5)  y2(r");  6)  y^Cm);  7)  z2(m);  8)  z^Cm).  °  1 


Table  2.  Initial  Condition  for  Trajectory  and  Filtering. 
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H 
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Key:  1)  dynamic  state;  2)  initial  condition  of  dynamic  state 
functions;  3)  initial  value  of  filtering;  4)  unit;  5-7)  m/sec; 
8-10)  1/sec;  11)  m. 


Table  3.  Standard  Deviation  for  Measurement  and  Process  Noise. 
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Key:  12)  symbol;  13)  numerical  value;  14)  noise  definition; 

15)  remark;  16)  0.0942/sec;  17)  2.0  m/sec  ;  18)  0.5  1/sec  ; 

19)  acceleration  measurement;  20)  angular  rate  measurement; 

21)  linear  acceleration  process;  22)  transverse  angular 
acceleration  process;  23)  15?  of  the  entire  processed  quantity. 


Table  4.  True  Values,  Filtering  Initial  Values  and  Filtering 
Results  for  Aerodynamic  Coefficients. 
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Key:  24)  aerodynamic  coefficient;  25)  real  value;  26)  initial 

values  of  filtering;  27)  results  of  filtering;  28)  optimum 
estimation. 


From  equations  (5)  and  (6),  and  the  relevant  data  in  Tables  /20 

1-4,  we  obtained  the  discrete  values  of  acceleration  and  angular 

rate  n'  ,  n„,  n  ,  p,q,r.  The  sampling  interval  was  chosen  to  be 
x  y  2 

0.02  second  as  shown  in  Figure  1.  These  data  are  stored  in 
the  machine  to  be  used  as  the  observation  data  for  filtering. 

Using  equation  set  (7)  and  the  relevant  data  in  Tables  1-4, 
the  optimum  estimated  value  cK  and  filtering  computation  square 
deviation  P^(+)  of  various  aerodynamic  coefficients  were  obtained 
based  on  filtering  iteration. 
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Figure  1.  Measurement  data. 

Key:  1)  p(l/sec);  2)  q  (1/sec);  3)  r(l/sec);  4)  t(sec) 


The  computed  results  are  shown  by  two  indicators;  one 
is  the  actual  estimated  error  Eg  of  filtering  and  the  second 
is  the  computed  error  of  filtering  Ef  .  The  formulas  are 
as  follows: 


(£«)*■*  ioo(c*~  C)/C 
(£,)«-  ±  ioov/P«(  +  )/1C| 


A 

where  C  „  represents  the  optimum  estimated  value  of  an  aero- 
dynamic  coefficient  at  the  K  observation  point,  C  represents 
the  actual  value  corresponding  to  the  aerodynamic  coefficient 
(see  Table  4),  and  PK(+)  represents  the  square  deviation  cor¬ 
responding  to  the  aerodynamic  coefficient  which  is  calculated 
from  the  filtering  equation  (7).  Basic  filtering  theory 
points  out  that  if  the  two  indicators  are  in  agreement  statis¬ 
tically  then  it  indicates  the  filtering  is  normal.  The  (E  )„ 
obtained  after  converging  can  be  used  as  the  actual  estimated 
error. 

The  last  3  columns  in  Table  4  give  the  filtering  results.  /21 

The  first  column  shows  the  optimum  estimated  values  when  K=50. 

The  two  latter  columns  represent  the  Eg  and  E^.  values  of  the 
13  aerodynamic  coefficients  corresponding  to  K=50  and  K=80, 
respectively. 

Figures  2-4  plot  the  filtering  process  of  7  major  aero¬ 
dynamic  coefficients .  The  solid  lines  in  the  figures  represent 
the  actual  estimated  error  E&.  The  symmetrical  dotted  lines 
are  the  filtering  calculated  error  Ef. 

Figure  5  shows  the  preliminary  extraction  results  of  a  certain 
actual  measured  data  which  corresponds  to  the  trajectory  re-en'ry 
condition.  In  the  figure,  the  analytical  solution  extraction 

rgi 

results  of  K>30  with  trajectory  measured  data  are  also  plotted.1-  - 

The  above  figures  reflected  the  following  facts: 

(1)  For  all  the  aerodynamic  coefficients,  regardless  of  whether 
it  is  the  actual  estimated  error  Eg  or  the  filtering  calculated 
error  Ef,  there  is  no  divergence  effect.  \E./E,\i&2 
For  the  major  aerodynamic  coefficients,  such  as 

CmC„,C*,  C„,  C.#,  CMl  ,  the  values  of  Eg  and  Ef  converge  very 
rapidly.  This  indicates  that  filtering  is  normal. 
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(2)  Through  a  50  point  filtering  computation  (which  is 
the  flight  data  in  1  second),  the  actual  estimated  error 
of  the  above  described  7  major  aerodynamic  coefficients 
converges  to  within  1.5%.  The  filtering  computation 
error  is  approximately  7%. 


-  E. 


Figure  2.  Filtering  Results  for  Aerodynamic  Force 
Coefficients.  ('  c.-  c,„  ). 

Key:  1)  t  (sec). 


Figure  3.  Filtering  Results  for  Aerodynamic  Moment 
Coefficients.  (Cma,  Cng). 


-  E. 


Figure  Filtering  Results  for  Control  Moment 
Coefficients  (cm6q,  cn6r). 

Key:  1)  t(sec);  2)  this  work;  3)  Reference  [81. 


For  coefficients  such  as  C„ rt  C,»,  C,tff  Cftr  etc.  , 

although  Ef  varies  convergently ,  yet  they  do  not  converge  as 
quickly  as  the  7  major  aerodynamic  coefficients  described 
above.  The  actual  estimated  error  Eg  does  not  have  any  sig¬ 
nificant  improvement.  This  indicates  that  under  given  conditions 
these  coefficients  do  not  significantly  affect  the  motion 
in  this  time  period. 

(3)  Preliminary  extraction  of  real  data  showed  that 
the  mathematical  model  and  filtering  program  in  this  paper 


could  adequately  suit  real  flight  data.  The  results  of  this 
work  are  very  close  to  those  of  the  analytical  solution  with 
trajectory  measured  data. 

In  addition,  several  points  are  explained  as  follows: 

(1)  Rigorously  speaking,  it  should  provide  many  initial 
values  of  xQ  randomly  and  through  Monte  Carlo  simulation, the 
final  statistical  average  value  is  given.  However,  the  above 
set  of  classical  situation  computation  is  already  capable  of 
reflecting  the  major  characteristics  of  filtering. 

(2)  For  the  real  flight  data,  to  provide  the  initial  square 
deviation  is  undoubtedly  an  important  problem.  The  ?q 
corresponding  to  an  aerodynamic  coefficient  can  be  given  based 

on  wind  tunnel  experiment  data  or  from  approximations.  The  PQ 

corresconaing  to  a  dynamic  state  variable  was  given  by  an 

f  81 

existing  extraction  method.  J 

(3)  Although  the  analysis  in  this  paper  has  the  various 
assumptions  described  in  the  previous  section,  this  basic 
method  may  be  extended  to  a  more  generalized  situation  such 
as  the  aerodynamic  non-linear  effect ,  non-exponent  type  of 
atmospheric  density;  various  types  of  measurement  errors, 
etc . 

4.  CONCLUSIONS 

1.  This  paper  analyzed  the  use  of  an  extended  Kalman 
filtering  method  to  determine  the  aerodynamic  coefficients 
of  a  re-entry  body  from  angular  rate  and  acceleration  data 
on-board  in  the  absence  of  trajectory  observation  data. 

Numerical  simulation  and  preliminary  real  data  extraction 
showed  that  the  solution  of  this  problem  is  practical. 


JKWIHMJ 


2.  As  long  as  no  stringent  initial  state  estimation  is 
required,  after  the  treatment  of  1-2  second  flight  data,  the 
error  of  a  major  aerodynamic  coefficient  can  converge  to 

a  small  range. 

3.  In  the  selection  of  parameters  and  analysis  of 
error,  more  in  depth  work  is  needed. 


We  wish  to  thank  Vice  Chief  Zhang  Hauxiu  and  Associate 

Professor  Zhang  Jinghuai  for  their  valuable  opinions. 
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DETERMINATION  OF  AERODYNAMIC  COEFFICIENTS 
FOR  A  RE-ENTRY  BODY  BY  MEANS  OF  AN 
EXTENDED  KALMAN  FILTER 


Jiang  Quanvoei,  Xu  Jinxhi,  Zhou  Shuying 
( Aerod ynanti c  Retearch  Center  of  Chino ) 

Abctract 

An  extended  Kalman  filter  is  used  to  determine  the  major  aerodynamic  co¬ 
efficients  of  a  re-entry  body  in  this  paper.  The  emphasis  is  pat  tm  estimating 
the  states  and  parameters  only  on  the  basis  of  the  re-entry  on-board  measure¬ 
ment  data  in  the  absence  of  trajectory  observation  data.  The  measured  data 
include  angular  rates  and  acceleration  obtained  from  the  rate  gyros  and  acce¬ 
lerometers  installed  in  the  vehicle.  A  mathematical  model  presented  is  based  u- 
pon  the  6-degree-of-freedom  motion  equations.  Both  ballistic  and  maneuvering 
re-entries  are  considered.  The  numerical  simulation  and  real  data  extraction 
sholr  that  the  presented  method  can  provide  satisfactory  results. 


A  PROGRAM  SYSTEM  FOR  DYNAMIC  ANALYSIS  OF 
AERONAUTICAL  STRUCTURES  -  HAJIF-II 
TEAM  FOR  DEVELOPING  HAJIF-II 

Institute  of  Aeronautical  Research  of  China 
Written  by  Guan  De  * 


ABSTRACT 

HAJIF-II  is  the  dynamic  analysis  system  for  aeronautical 
structures  developed  under  the  leadership  of  the  Institute  of 
Aeronautic  Research  of  China.  It  is  capable  of  carrying  out 
the  calculation  of  natural  dynamic  characteristics  of  structures 
as  well  as  the  calculation  of  aircraft  flutter  with  an  active 
control  system  and  the  calculation  of  gust  response.  It  has 
31  fixed  flow  routes  and  2600  FORTRAN  statements.  It  allows 
the  use  of  99  elementary  substructures  and  each  substructure 
can  have  7000  degrees  of  freedom.  In  the  flutter  and  gust 
response  computation,  it  is  possible  to  use  50  modes.  The 
panel  number  in  nonsteady  aerodynamic  calculation  can  reach 
300.  In  the  management  of  the  stiffness  matrix  and  mass  matrix, 
the  hypermatrix  method  is  used  for  "macroscopic  treatment" 
and  the  effective  column  method  is  used  for  "microscopic 
treatment"  to  develop  a  new  simultaneous  iteration  algorithm 
to  improve  the  efficiency  of  the  calculation  of  real  characteristic 
values.  A  more  complete  state  synthesis  method  computation 
program  has  been  designed.  In  addition  to  fixed  and  free 
interface  methods,  multilevel  synthesis  and  step-by-step  com¬ 
putation,  a  curve-fitting  method  is  introduced  to  transform 
the  harmonic  oscillation  aerodynamic  force  into  the  Laplace 
plane  technique.  This  system  has  been  used  to  calculate 
a  number  of  typical  aircraft  structures.  Good  results  are 
obtained. 

•Received  on  December  15,  1981. 
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1.  INTRODUCTION 

KAJIF-II  is  the  aeronautical  structure  dynamic  analysis 
system  developed  by  the  Institute  of  Aeronautics  Research  of 
China.  Its  function  includes  the  calculation  of  the  natural 
dynamic  characteristics  of  structures,  the  calculation  of  aircraft 
flutter  with  an  active  control  system  and  the  calculation  of 
gust  response. 

In  the  development  of  the  system,  a  series  of  measures 
has  been  adopted  to  satisfy  generality,  flexibility,  reliability, 
automation,  high  efficiency,  expandability,  ease  of  correction, 
and  diagnostic  capability  requirements.  Furthermore,  special 
attention  has  been  paid  to  the  guiding  principle  of  advancing 
on  the  basis  of  practicality. 

2.  CAPABILITIES  AND  SIZE 

In  the  calculation  of  structural  natural  dynamic  characteristics 
of  HAJIF-II,  it  is  allowed  to  use  99  elementary  substructures  and 
each  substructure  has  no  more  than  7000  degrees-of-freedom.  In 
nonsteady  aerodynamic  calculation,  the  total  number  of  panels 
can  reach  300.  In  the  flutter  and  gust  response  calculations, 
it  is  possible  to  use  50  modes. 

HAJIF-II  provides  31  fixed  flow  routes.  That  is,  there  are 
7  computational  flows  in  the  calculation  of  the  natural  dynamic 
characteristics  of  structures  (whole  structure,  fixed  interface 
method  single  synthesis,  free  interface  method  single  synthesis, 
fixed  interface  method  multi-level  synthesis,  free  interface 
method  multi-level  synthesis,  fixed  interface  method  step-by- 
step  synthesis,  free  interface  method  step-by-step  synthesis), 

2  flows  in  the  flutter  calculation  flow ,  (v-g  method,  p-k  method), 

1  gust  response  calculation  flow,  and  21  combination  flow  rates  to 
calculate  natural  dynamic  characteristics ,. flutter ,  and  gust  response. 
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3-  STIFFNESS  AND  MASS  MATRIX 


In  the  natural  dynamic  characteristics  calculation  of  / 2 6 

structures  of  HAJIF-II,  a  finite  element  displacement  method 
is  used.  A  structure  can  be  discretized  into  a  finite  element 
model,  or  can  be  approximated  by  a  single  beam.  A  structure 
can  be  analyzed  as  an  entity,  or  can  use  the  mode  synthesis 
method. 


In  the  element  warehouse  of  HAJIF-II,  there  are  11  tyres 
of  bar,  plate,  beam,  and  film  elements  to  simulate  the  wir.g, 
body,  external  attachment,  and  their  correct ior.s  . 

In  order  to  improve  the  efficiency  of  calculation,  the 
following  measures  are  taken: 

(1)  The  technique  to  vary  the  degree  of  freedom.  Regardless 
of  whether  it  is  an  elementary  or  higher  level  substructure, 

the  program  automatically  determines  the  real  degree  of  freedom 
of  each  nodal  point  to  eliminate  ineffective  degrees  of 
freedom. 

(2)  The  Cuthill-Meckec  method^-  is  used  in  the  re-labeling 
of  multi-level  substructural  nodal  points  and  the  optimized 
utilization  of  band  width.  In  addition,  some  special  treatments 
are  also  performed. 

(3)  A  compact  assembly  method  is  used  to  only  assemble 
non-zero  nodal  point  panels. 

(4)  The  improved  hypermatrix  technique.  In  an  usual 
hypermatrix  method,  there  are  considerable  zero  elements  in 
a  non-zero  submatrix.  For  this,  HAJIF-II  first  conducts  a 
"macroscopic  treatment"  to  the  submatrix  in  the  hypermatrix 
method,  and  then  carries  out  a  "microscopic  treatment".  For 
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the  stiffness  matrix,  the  effective  column  concept  is  used  to 
eliminate  the  useless  zero  elements  in  a  non-zero  submatrix. 
For  a  mass  matrix,  because  it  is  unchanged  in  the  entire 
iteration  process,  a  nodal  point  control  method  is  used. 

The  improved  hypermatrix  method  can  effectively  reduce  the 
storage  requirement  as  shown  in  Table  1. 


Table  1.  Comparison  between  the  internal  storage  used  by  HAJIF-II 
and  that  by  the  conventional  hypermatrix  technique. 


1 — 

*  n 

■a — 

»  ft 

i 

if  «««**£* 

S'  ****#£> 

3  -mmamm  J 

HAJIF-I 

HAJ1F- 1 

i 

170 

107712  I 

70SS7 

21312 

1504 

i 

210 

11520  1 

0500 

0500 

1030 

s 

130 

0320  j 

1450 

5357 

070 

Key:  1)  example;  2)  order;  3)  conventional  hypermatrix; 

4)  storage  of  the  stiffness  matrix;  5)  storage  of  the 
mass  matrix;  6)  conventional  hypermatrix. 


In  the  constraint  treatment  area,  HAJIF-II  can  perform 

single  point  or  multiple  constraint  treatment  with  regard 

f2l 

to  the  degree  of  freedom  of  a  nodal  point L  J.  The  function 
of  the  former  is  to  designate  the  displacement  values  of  certain 
degrees  of  freedom  to  simulate  the  symmetric  and  asymmetric 
conditions  and  boundary  conditions,  as  well  as  to  conduct  zero 
stiffness  direction  treatment.  The  latter’s  function  is  to 
express  the  displacements  of  part  of  the  degrees  of  freedom  as 
the  linear  combination  of  the  displacement  of  other  degrees 
of  freedom  in  order  to  treat  the  coordinating  and  stiffness 
elements  between  the  degrees  of  freedom. 
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IV.  ANALYSIS  OF  HEAL  CHARACTERISTI C  VALUES 

In  the  natural  dynamic  characteristics  calculation  of  structures, 
the  following  equation  is  used 

An-Xmu  (1) 

where  k  and  m  are  the  stiffness  matrix  and  mass  matrix,  respectively, 

X  and  u  are  characteristic  value  (circular  frequency  square)  and 
characteristic  vector  (natural  mode),  respectively. 

In  order  to  improve  the  calculation  efficiency,  a  new 
algorithm  is  developed  on  the  basis  of  theoretical  analysis  of 
existing  simultaneous  iteration  methods  and  large  amounts  of 
numerical  tests  by  combining  the  advantages  of  these  methods. 

Its  special  point  is  that  the  mass  matrix  is  not  decomposed. 

Rather,  a  "dimension  lowering"  treatment  and  other  improvements  /27 

are  used.  Therefore,  compared  to  the  existing  algorithms, 
this  new  algorithm  has  the  same  converging  rate.  However,  the 
calculation  load  in  each  iteration  step  is  the  least  and  the 
efficiency  is  high  as  shown  in  Table  2. 

When  using  projection  matrix  to  solve  the  characteristic 

value  problem,  the  QR  or  Jacobi  method  is  used.  The  Chebyshev 

polynomlnal  or  origin  displacement  is  used  to  perform  acceleration 

r  q  1 

treatment.  The  characteristics  of  the  Sfurm  seriesLJJ  are 
used  to  carryout  the  missing  root  check.  This  can  be  carried 
out  in  combination  with  the  determination  of  the  number  of 
trial  vectors.  The  wave  array  element  elimination  method  with  a 
buffer  is  used  to  solve  the  matrix. 
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Table  2.  Comparison  of  efficiencies  of  various  simultaneous 
iteration  algorithms. 
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Key:  1)  example;  2)  order;  3)  half  band  width;  4)  calculation 
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V.  .MODE  SYNTHESIS 

In  HAJIF-II,  a  more  complete  mode  synthesis  program  has 
been  designed.  In  addition  to  the  two  major  methods  —  fixed 

fil  El 

interface  and  free  interface  methodL  ’  ,  we  also  developed 

a  level  synthesis  and  step-by-step  synthesis  method  of  our 
own  to  form  six  types  of  synthesis  methods  as  described  in  Section 

2. 


In  the  free  interface  method,  we  should  consider  the 
residual  stiffness  and  residual  inertia  effects  in  higher 
order  modes.  Through  trial  calculation,  it  was  discovered 
that  the  residual  inertia  effect  could  be  neglected.  Hence, 
in  HAJIF-II  only  the  algorithm  considering  the  residual  stiff¬ 
ness  effect  is  used. 

The  free  interface  method  of  HAJIF-II  processes  the 
attached  mode  into  a  "quasi  constraint  mode"  and  the  free 

r  c  I 

interface  main  mode  into  a  "quasi  fixed  interface  main  mode"1  J 
so  that  the  calculation  format  of  the  free  interface  method 
is  unified  with  that  of  the  fixed  interface  method. 
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In  the  search  for  the  softness  matrix  of  the  free  interface 
main  mode  and  free  structure,  the  "imaginary  structure"  concept 
is  used.  The  stiffness  matrix  of  the  imaginary  structure  is 
k*S£k+m,  and  the  mass  matrix  is  still  m.  Hence,  the  softness 
matrix  of  a  free  structure  is  (k*)”^. 

The  differences  among  single  synthesis,  multi-level 
synthesis  and  step-by-step  synthesis  are  shown  in  Figure  1. 

The  multi-level  synthesis  is,  in  principle,  an  extension  of  the 
single  synthesis.  However,  in  each  level  of  synthesis  in  the 
step-by-step  synthesis  method,  tnere  are  two  substructures 
participating.  Among  them,  one  substructure  participates 
using  the  reduced  generalized  degree  of  freedom  and  the  other 
substructure  directly  participates  with  the  physical  degree 
of  freedom. 

In  order  to  suit  the 
requirements  of  multi-level 
synthesis  and  step-by-step 
synthesis,  some  development  is  needed 
with  regard  to  frequency  selec¬ 
tion  judgement  in  HAJIF-II. 

HAJIF-II  can  treat  rigid 
substructures.  Only  when 
using  multi-level  synthesis, 
it  is  required  to  process  it 
during  the  first  synthesis. 

Figure  1.  Models  for  cal¬ 
culation  of  various  modal 
synthesis  techniques. 

Key:  1)  single  synthesis; 

2)  multi-level  synthesis; 

3)  step-by-step  synthesis 
(synthesized  using  the  physi¬ 
cal  equation  for  those  with¬ 
out  the  a  symbol  ) . 


VI.  NONSTEADY  AERODYNAMIC  FORCE  CALCULATION 
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Serious  flutter  of  modern  aircraft  mostly  occurs  in  the  near 
sonic  speed  region  at  low  altitude.  For  wing  surface  with  a 
medium  or  small  aspect  ratio,  the  results  calculated  based  on 
the  subsonic  nonsteady  aerodynamic  theory  could  better  reflect 
the  basic  flutter  characteristics.  In  HAJIF-II  we  first  esta¬ 
blished  a  subsonic  horseshoe  vortex  —  oscillating  doublet 
r  7 1 

mesh  method1- '  J  to  treat  multiwings  in  space. 

In  order  to  calculate  the  non-steady  aerodynamic  force, 
it  is  necessary  to  insert  the  vibrational  mode  at  the  structural 
nodal  point  to  the  ooint  of  aerodynamic  force.  HAJIF-II  uses 

r  8 1 

the  two-way  spline  curves  function  method1  . 

r  q  i 

In  HAJIF-II,  a  curve  fitting  method  used  by  Rogai  was 
Introduced  which  utilizes  the  resonating  aerodynamic  data 
to  find  the  approximate  expression  of  non-steady  aerodynamic 
force  on  the  Laplace  plane. 


VII.  FLUTTER  CALCULATION 


In  HAJIF-II,  v-g  and  p-k  methods  are  used  to  carryout 
flutter  calculation^10-^.  The  v-g  method  uses  the  following 
flutter  equation 


(2) 


The  p-k  method  uses  the  following  flutter  equation  to  be 
taken  Into  account  by  the  active  control  system 
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where  "K,  H,  and  A  are  generalized  stiffness,  mass,  and  aero¬ 
dynamic  force,  respectively;  q,a,  6  are  the  generalized  coordinates 
of  the  HAJIF-II  users  to  add  and  control  modes;  p  and  V  are  the 
atmospheric  density  and  flight  velocity;  X  is  the  complex 
characteristic  value;  bg  is  a  reference  length. 

HAJIF-II  uses  LR,  QR,  and  reverse  root  method  with  origin 
displacement  to  solve  the  complex  characteristic  value  problem. 

They  are  available  for  the  users  to  choose. 


VIII.  gust  response  calculation 


KAJIF-II  uses  the  following  equation  to  calculate  the 
aerodynamic  elastic  frequency  response^] 


+  Ca«(  Oa.#<  *  »  ){  0  }--J  pr*U,(0> 


<*0 


tfn-£5fry<'» 


where  R/D  is  the  transfer  function  of  the  active  control  system. 
The  footnote  g  represents  gust  and  s  is  the 


In  the  calculation  of  gust  response,  it  is  possible  to 
use  the  Karman  power  spectrum  or  to  provide  the  users  with  a 
fixed  form. 


LX.  PROGRAM  ORGANIZATION  AND  USER  INTERFACE 

HAJIF-II  uses  a  two  level  management  system.  The  first 
level  is  the  monitoring  control  program  which  is  composed  of 
a  series  of  dispatching  languages.  The  second  level  is  the 
processing  program.  In  fact,  it  is  dispatched  by  the  monitoring 
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control  program,  which  is  the  highest  level  subroutine  capable 
of  finishing  a  certain  mathematical  or  mechanical  function  inde¬ 
pendently.  It  is  not  possible  to  be  mutually  dispatched  between 
two  processing  programs.  The  input  and  output  of  each  processing 
program  are  only  related  to  data  documentation. 

In  the  monitoring  control  program  of  HAJIF-II,  a  common 
number  group  with  a  variable  length  is  established  to 
be  given  to  various  processing  programs.  In  addition, to  be  able 
to  establish  a  very  small  amount  of  the  number  group  by  the 
processing  program,  other  numbers  are  formed  completely  from 
the  above  mentioned  common  number  group. 

For  a  calculation  requiring  a  larger  internal  storage,  the 
system  can  automatically  judge  whether  external  storage  is 
necessary. 

The  entire  data  of  HAJIF-II  is  called  the  "dictionary". 

The  dictionary  is  formed  by  several  volumes  and  the  number 
of  volumes  can  be  expanded.  Each  volume  has  100  books  and 
each  book  is  made  up  of  10 0  chapters.  The  volume,  in  principle, 
is  divided  according  to  the  processing  program.  The  major 
purpose  of  a  book  Is  to  facilitate  the  performance  of  some 
repeated  calculation.  A  chapter  is  the  entire  data  of  a 
"read"  or  "write"  statement,  which  is  the  basis  of  data  manage¬ 
ment  . 

The  document  management  system  of  HAJIF-II  can  automatically 
conduct  document  distribution  and  correlated  variable  control. 

As  long  as  the  input,  output  table  and  the  chapter  address  are 
provided,  the  system  can  independently  complete  reading  and 
writing  work. 
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In  HAJIF-II,  user- oriented  broad  command  languages  were 
designed.  One  is  used  to  initiate  leading  fixed  flow  route 
and  the  other  one  is  used  to  dispatch  21  processing  programs 
for  the  users  to  form  a  "user  flow".  In  view  of  the  fact  that 
flow  organization  of  mode  synthesis  is  more  complicated,  we 
also  designed  the  special  composite  language. 

In  order  to  facilitate  the  preparation  of  data,  the 
data  generation  system  of  HAJIF-II  allows  the  use  of  two 
input  forms:  one  is  numerical  data  input  and  the  other  is 
a  mixed  input  of  numerical  values,  data,  and  topological 
description. 

Y 

The  function  of  topological  description  includes  the 
simple  description  of  nodal  point  coordinates  of  ideal  parts 
and  element  information,  data  abbreviation  with  regularity, 
data  correction,  addition  and  elimination,  etc. 

As  for  the  error  diagnosis  of  HAJIF-II,  the  focus  is 
placed  on  the  checking  of  the  original  data.  It  is  capable 
of  providing  over  one  hundred  types  of  information  according 
to  the  requirement  of  the  users. 

X.  EXAMPLE 

during  the  development  of  HAJIF-II,  over  ten  examples  were 
calculated.  Now,  let  us  introduce  three  types  of  combined 
examples . 

(1)  Calculation  <_  the  Natural  Dynamic  Characteristics 
and  Flutter  of  a  Triangular  Wing.  The  wing  is  simulated 
by  bars  and  shearing  plates.  The  fuselage  and  flat  tail  are 
simulated  by  beams,  and  the  wing- fuselage  connection  is  simulated 
by  transition  beams.  There  are  315  nodal  points  and  1011  degrees 
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of  freedom.  The  comparison  of  calculated  natural  dynamic  charac 
teristics  and  the  resonance  test  results  is  shown  in  Figure  2. 

It  coincides  very  well.  In  the  meantime,  v-g  and  p-k  methods 
are  used  to  calculate  flutter  and  the  results  agree  with  the 
wind  tunnel  test  results. 

(2)  The  calculation  of  modal  synthesis  of  the  flat  wing  — 
fuselage  —  attachment  combination.  The  structure  is  simulated 
entirely  by  beams.  There  are  nodal  points,  138  degrees  of 
freedom,  and  4  substructures  as  shown  in  Figure  1.  The 
calculated  natural  frequencies  are  shown  in  Table  3.  It  is 
obvious  that  all  the  synthesis  methods  have  excellent  accuracy. 

(3)  Calculation  of  Gust  Response  of  B-47  Airplane.  The 
computation  data  in  Reference  [11]  was  used.  The  wing  and 
fuselage  are  simulated  by  beams  and  the  flat  tail  is  a  rigid 
body.  In  the  aerodynamic  force  calculation,  the  interference 
between  the  wing  and  the  flat  tail  was  included.  In  the 
calculation  of  the  response,  the  atmospheric  turbulence  power 
spectrum  measured  in  flight  tests, provided  in  Reference  [11],  was 
used.  The  results  are  shown  in  Figure  3*  Because  the  nonsteady 
aerodynamic  force  used  in  Reference  [11]  is  different  from  the 
one  used  in  HAJIF-II,  the  results  of  these  two  calculations 

are  not  quite  the  same. 
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Figure  2.  Comparison  of  calculated  natural  frequencies 
and  modes  with  resonance  test  results. 

Key:  1)  mode;  2)  first  level  wing;  3)  second  level  wing; 
4)  third  level  wing;  5)  nodal  line;  6)  experimental; 

7)  calculated;  8)  calculated  frequency;  9)  experimental 
frequency. 


able  3.  Comparison  of  various  synthesis  techniques. 
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Key:  10)  mode;  11)  whole  structure;  12)  fixed  interface  single 

synthesis;  13)  free  interface  single  synthesis;  1M)  fixed  interface 
multilevel  synthesis;  15)  free  interface  multilevel  synthesis; 

16)  fixed  interface,  step-by-step  synthesis;  17)  free  interface, 
step-by-step  synthesis. 
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Figure  3.  Comparison  of  calculated  gust  responses  with 
the  results  of  Reference  [11]. 

Key:  1)  acceleration  response  amplitude  at  center-of- 

gravity  ( foot /sec  )/( foot /sec ) ;  2)  frequency  (Hertz); 

3)  rigid  axis;  4)  Reference  [11]. 


CONCLUSION  REMARKS 

This  paper  briefly  described  the  method  and  function  of 
the  HAJIF-II  dynamic  analysis  system  for  aeronautical  structures. 
The  system  incorporated  several  new  measures  and  quite  good 
trial  calculation  results  were  obtained.  The  system  will  be 
further  developed.  Due  to  the  limitation  in  space,  this  paper 
might  appear  to  be  too  oversimplified.  The  authors  wish  to 
apologize  to  the  readers  in  this  regard. 
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A  PROGRAM  SYSTEM  FOR  DYNAMIC  ANALYSIS  OF 
AERONAUTICAL  STRUCTURES 
(HAJIF-II) 

Team  for  Developing  HAJIF-ll 

Penned  by  Guan  De 
Abstract 

HAJIF-ll  is  a  program  system  developed  by  Chinese  Aeronautical  Estab- 
lishment(CAE).  It  is  able  to  accomplish  the  calculation  of  the  modal  charac¬ 
teristics  of  aircraft  structures  as  well  as  the  flutter  and  gust  response  analysis 
with  the  active  control  system  taken  into  account.  T|se  structural  model  may 
be  composed  of  99  substructures  each  with  70QO  degrees  of  freedom.  300 
panels  may  be  used  in  the  calculation  of  nonsteady  aerodynamic  forces  and  50 
modes  in  the  flutter  and  gust  response  analysis.  The  data  generation  system 
permits  the  flexible  use  of  numerical  data  and  topological  description.  31  pre¬ 
scribed  computational  flows  are  supplied  and  an  user  can  also  organize  his  own 
computational  flows  as  he  needs.  A  structure  can  be  discretized  into  finite 
elements  or  simulated  by  single  spar.  For  the  management  of  the  stiffness  and 
mass  matrix  a  modified  hypermatrif  method  is  employed  to  omit  all  of  inacti¬ 
ve  zpro  elements  more  effectively.  A  new  algorithm*  called  a  revised  simulta¬ 
neous  iteration  procedure,  has  been  developed  to  solve  the  real  eigenvalue  pro¬ 
blem  and  is  more  effective  than  the  current  agorithm.  Modal  synthesis  techni¬ 
que  with  both  free  and  fixed  interfaces  is  adopted.  Besides,  two  new  methods 
of  synthesis  have  been  developed  from  the  concept  of  multilevel  substructures. 
Nonsteady  aerodynamic  forces  are  calculated  by  means  of  subsonic  doublet  latti¬ 
ce  method  for  multiwings  and  aerodynamie  forces  in  Laplace  plane  can  be  ap¬ 
proximated  'with  a  curve-fitting  procedure  based  on  sinusoidal  data.  Flutter 
equations  are  solved  by  V'-g  and  p~k  methods  and  the  continuous  atmosphe¬ 
re  turbulence  are  used  in  the  gust  response  analysis.  The  system  consists  of  a 
sequence  of  functional  modules  so  it  can  be  modified  and  extended  easily.  An 
advanced  file  management  system  has  also  been  developed.  There  are  app¬ 
roximately  26000  FORTRAN  IV  statements  in  the  system-  The  HAJIF-II  was 
applied  to  analyzing  a  number  of  typical  aircraft  structures  and  gave  good 
results. 


SYSTEM  IDENTIFICATION  AND  AIRCRAFT  FLUTTER 
Zhang  Lingmi * 

Nanjing  Aeronautical  Institute 


ABSTRACT 

After  briefly  describing  the  application  of  system 
identification  technique  to  the  study  of  aircraft  flut¬ 
ter,  this  paper  introduced  three  types  of  system  identi¬ 
fication  methods.  The  emphasis  was  placed  on  the  essence 
of  the  algorithm  and  its  theoretical  basis.  The  relevant 
data  processing  problem  was  also  properly  explained.  The 
transfer  function  method, based  on  complex  modal  analysis 
and  the  optimization  technique, could  accurately  identify 
all  the  modal  parameters.  It  was  successfully  used  in 
the  flutter  model  experiment  of  a  wing  with  external  bodies. 
This  method  could  be  extended  to  the  conditions  for  wind 
tunnel  and  flight  tests  with  response  data  only.  The  cor¬ 
relation  and  least  square  techniques  were  adopted  in  the 
impulse  response  and  auto-regressive  moving  average  model 
method  in  order  to  obtain  accurate  results  under  strong 
interference.  The  latter  could  also  obtain  the  mean  square 
deviation  between  the  estimated  value  and  the  actual  value. 
Both  time  domain  methods  could  directly  use  measured  samp¬ 
ling  data  to  perform  system  identification  without  special 
signal  analysis  equipment. 


I .  INTRODUCTION 

The  so-called  system  identification  in  structural  dynamics  esta¬ 
blishes  a  mathematical  dynamic  model  according  to  the  response  (out¬ 
put)  of  a  known  excitation  (input)  for  a  structural  system,  including 
the  modal  parametric  model  in  the  generalized  coordinate  system  and 
the  structural  parameters  model  In  the  physical  coordinate  system. 

The  more  practical  method  is  to  provide  the  input,  output  data  of  an 
actual  system  to  determine  a  model  structure  and  Its  parameters, 
using  a  certain  optimization  method  based  on  a  type  of  model  so  that 
they  are  best  matched  with  the  system  under  certain  guide  lines.  For 
a  discretized  linear  steady  system,  the  so-called  model  structure  con¬ 
sists  of  the  determination  of  the  system  (or  degree  of  freedom).  Once 
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the  order  is  determined,  system  identification  is  transformed  into 
parameter  identification  or  parameter  estimation  problems.  In  a 
complete  system  identification  process,  in  addition  to  the  structure 
determination  and  parameter  estimation  parts,  it  is  also  necessary 
to  carry  out  model  verification  (for  a  linear  system  the  main  item 
is  the  verification  of  order)  and  parametric  mean  square  deviation 
analysis.  In  order  to  further  improve  the  effectiveness  of  identifi¬ 
cation,  experimental  design  must  also  become  one  aspect  in  the  theor¬ 
etical  study  on  identification  [1]. 

In  the  recent  decade,  the  ideas  and  methods  of  system  identifi¬ 
cation  have  been  fairly  successful  in  structural  dynamics.  It  has 
promoted  the  ground,  wind  tunnel  and  flight  tests  of  aircraft  signi¬ 
ficantly.  In  the  ground  and  wind  tunnel  tests  of  models  and  flight 
tests  of  aircraft,  it  is  possible  to  more  accurately  identify  the  rel¬ 
evant  modal  parameters  [2,3].  For  a  windfall  flutter,  the  method  rely¬ 
ing  on  the  extrapolation  of  damping  to  determine  the  critical  point 
has  been  proven  ineffective.  However,  the  adaptation  of  the  system 
identification  method  can  determine  the  coefficients  of  the  flutter 
differential  equation  based  on  experimental  data.  From  these  coeffi¬ 
cients,  it  is  possible  to  more  accurately  obtain  the  flutter  critical 
velocity  [4,5].  For  the  ground  vibration  tests  of  the  entire  plane, 
the  determination  of  the  proper  adjustable  force  is  very  difficult. 
Presently,  the  major  solution  is  to  identify  the  transfer  function 
matrix  based  on  the  experimental  data  of  a  single  point  excitation 
and  the  shaker  adjustment  force  of  the  pure  mode  [6,7].  In  the  mean¬ 
time,  the  techniques  which  use  the  single  point  transient  or  random  exci¬ 
tation  transfer  function  to  identify  the  modal  parameters  of  the  whole 
aircraft  have  been  successfully  developed  [6,8]. 

System  identification  techniques  not  only  can  be  used  in  the 
ground,  wind  tunnel  and  flight  vibration  tests,  but  also  can  further 
serve  as  flutter  analysis  and  synthesis.  From  the  modal  parameters 
of  the  aircraft  (parts  or  the  whole  plane),  it  is  possible  to  further 
identify  the  stiffness  and  mass  matrices  of  structures  to  be  used  for 
structure  modification  and  model  optimization  [9,10].  For  an  aero¬ 
dynamic  elastic  model,  usually  only  numerical  solutions  can  be  given 
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at  the  present  time.  Using  the  identification  method,  it  is  possible 
to  obtain  the  mathematical  expression  in  a  rational  form  [11]  to 
facilitate  the  analysis  and  synthesis  of  the  flutter  of  the  stabi¬ 
lizer  system  with  control  and  the  flutter  active  control. 

This  paper  mainly  introduces  three  identification  methods  suit¬ 
able  for  studying  flutter,  i.e.,  the  transfer  function  method,  impulse 
response  method  and  the  auto-regressive  moving  average  (ARMA)  model 
method.  The  emphasis  is  placed  on  the  essence  of  the  identification 
algorithm  and  its  theoretical  basis.  With  regard  to  the  relevant 
data  processing  problem,  it  is  also  properly  explained.  These  methods 
are  not  only  suitable  for  the  sine  experimental  data,  but  also  trans¬ 
ient  and  random  experimental  data,  as  well  as  the  condition  of  using 
wide  frequency  band  random  natural  excitation  with  response  data 
available  alone. 

II.  TRANSFER  FUNCTION  METHOD 

After  the  discretization  of  a  linear  steady  system.  Its  dynamic 
equation  can  be  described  by  a  matrix  differential  equation  In  the 
following: 

CA/l<X>  +  CC)<i>+C/0<*>-{/(0>  (1) 

{/(Oh  {*(<)>  are  the  N-dimensional  external  force  and  the  displace¬ 
ment  response  vector,  respectively.  (A/),  CJChCCl  are  the  Nth  order 
mass,  stiffness  and  damping  matrices,  respectively.  [M]  is  orthogonal 
and  [K]  and  [C]  are  orthogonal  or  semi-orthogonal.  N  is  the  discret¬ 
ized  degree  of  freedom. 

By  conducting  Laplace  transformation  on  both  sides  of  the  equa¬ 
tion  and  assuming  the  initial  condition  is  zero,  we  obtain 

(**(Af!3  +  *  CC}  +  C/0)UI*)>-<F(«»  (2) 

where  (F(*)h  {A((«)>  are  the  Laplace  transform  of  {/(<))  and  '•XX  <  )> 
respectively;  s  is  the  Laplace  multiplier  (complex  number).  Then, 
equation  (2)  can  be  expressed  as 

(3) 

cz(*))ur(o>-{F(*)> 

Ctf(*)XF(  «)>-{*(«)> 

(4) 


or 
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where 


cz(OWtMi+*cc)+c/o  (5) 

is  the  generalized  dynamic  stiffness  matrix.  Its  inverse  matrix 

C Hi  *  ))-CZ(  *  )r'-adiCZ(  *  ))/detCZ(  *  ))  (6  ) 

is  the  generalized  dynamic  softness  matrix  or  the  transfer  function 
matrix. 

Let  sr  be  the  characteristic  root  corresponding  to  the  equation 
(3)  and  let  {1^}  be  the  corresponding  characteristic  vector,  i.e., 

CZO,)MtM-  0  (7) 

For  a  sub-critical  damping  system,  sr  and  {^r>  are  complex  numbers 
which  are  complex  conjugates  appearing  in  pairs. 


From  equation  (6),  we  know  that  the  transfer  function  can  be 
expressed  as  a  rational  function  of  s.  Furthermore,  it  can  be  expanded 
around  the  characteristic  root. 


tr 


(«<*»- 


(8) 


r  *  1 


where  [Arl  is  the  residual  number  matrix  corresponding  to  sr-  The 
symbol  *  represents  the  complex  conjugate.  It  can  be  proven  that  [Ar] 
and  the  characteristic  vector  { 4»r }  have  the  following  analytical 
relation  [12,13] 


P,-<^}TC2*rCA/I  +  (CD3{W 


(9) 


When  the  system  damp  is  the  structure  damp  or  proportionality 

* 

damp,  the  complex  vibrational  vector  {^r>  becomes  the  real  vibration¬ 


al  vector  {$>r) 


From  equation  (9),  we  can  get 
(  i  —V'—  1 ,  0r-/„(ir) 

(  m„  is  the  real  mode  mass  of  the  r-th  order).  A 


is  also  a  pure  imag¬ 


inary  number.  After  substituting  into  equation  (8),  we  obtain 

CH(Ol»  2  ■  Sn'+sC'+KT  ao) 

f  •  1 


Kr  and  Cr  are  the  r-th  order  real  mode  stiffness  and  damp,  respectively, 
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s  =:'  co 


By  letting  #,«■ -<x,+jpr,  j4,=ur+iv„  and  substituting  *  =:co  into 
the  general  relation  of  the  transfer  matrix  of  equation  (8),  we  can 
get  the  analytical  expression  between  the  frequency  repsonse  and  the 
modal  parameters 

(u.+i'u,}  ,  (a,— fry] 

(11) 


V  .  .  Ca,-iaO 

ttf  («<■>))  2^  ar+i(a>  —(},)  +  a,+  <(<»+0,) 


A  row  (corresponding  to  a  single  coordinate  excitation  multi-coordi¬ 
nate  measurement)  or  a  column  (corresponding  to  the  single  point  exci¬ 
tation  single  coordinate  measurement  of  each  respective  coordinate) 
of  C can  determine  the  whole  modal  parameters. 


Note  that  (ft)  is  the  measured  frequency  response  data  correspond¬ 
ing  to  various  frequencies  tested  and  (&)“(<*,,  0„  p, . aN,  0K,  uK,  »W)T 

is  the  modal  parameter  vector  to  be  identified.  Hence,  the  total  square 
deviation  /(  8  )«■■{//(  8  )  —  ff}TlH(8)  —  of  the  frequency  response  data 
and  model  value  {//( 8  )>  can  be  used  as  the  identification  guideline 
function.  Note  that  H(0)  is  the  nonlinear  function  of  the  parameters 
to  be  determined  (such  as  ar,  g^).  Therefore,  the  modal  parameter 
identification  problem  is  transformed  into  the  optimization  of  J(0)  = 
min  which  can  be  solved  by  iteration:  to  first 

give  the  initial  value  of  the  parameter  {0}  to  be  determined  and  to  use 
the  nonlinear  least  square  multiplication  method  to  obtain  the  incre¬ 
ment  { A 0 >  as  the  optimizing  direction  and  then  to  use  an  extrapolation 
method  to  find  the  step  length  factor.  The  actual  algorithm  can  be 
referred  to  in  [13]. 


The  row  data  of  parameter  identification  using  the  transfer  func¬ 
tion  method — frequency  response,-  can  be  obtained  by  two  ways:  one  is 
to  use  steady  state  sinusoidal  excitation,  tracking  filter  data  anal¬ 
ysis,  and  the  other  is  to  use  wide  frequency  band  excitation  (fast 
sinusoidal  scanning  frequency,  impulse,  pure  random,  pseudo-random  or 
periodical  random,  etc.)  fast  Fourier  transfer  (FFT)  data  processing. 
The  frequency  response  of  the  former  represents  the  complex  amplitude 
ratio  of  response  and  excitation  force  and  the  latter  represents  the 
ratio  of  the  mutual  power  spectrum  and  the  individual  power  spectrum 
of  input  and  output. 


The  main  advantage  of  steady  state  sinusoidal  excitation,  track¬ 
ing  filter  analysis  is  the  high  data  accuracy.  For  wind  tunnel  model 
flutter  tests  on  wings  with  external  bodies,  the  accuracy  of  the 
modal  frequencies  identified  could  reach  1% .  The  vibrational  model 
orthogonality  tests  could  reach  the  10?  level  required  by  the  finite 
element  mathematical  model  actually  used  in  the  analysis,  testing 
and  modification  of  flutter  [1*1]. 


The  major  advantages  of  the  wide  frequency  bandwidth  excitation, 
FFT  analysis  technique  are  the  fast  testing  speed  and  the  adaptation 
to  the  wind  tunnel  and  flight  test  conditions.  The  disadvantages  are 
the  small  signal-to-noise  ratio  and  the  low  frequency  resolution 
which  affect  the  data  accuracy.  Presently,  measures  to  improve  the 
accuracy  of  data  have  been  developed,  such  as  the  power  spectrum  data 
smoothing  and  bandwidth  selectable  FFT  analysis  (BSFA),  etc. 


In  the  subcritical  wind  tunnel  and  flutter  tests,  sometimes  it  is 
possible  to  directly  use  natural  bandwidth  random  excitation  (such  as 
turbulence).  At  this  time,  the  input  is  a  multiple  input  which  cannot 
be  directly  measured.  Assuming  that  the  external  force  frequency  spec¬ 
trum  F(w)  near  the  modal  frequency  is  a  constant,  then  the  response 
spectrum  can  be  expressed  asr  « 

Cu, +»'»,) 


{U'  +  iV,} 


N 

2  a, 4-  »  (ao  -pj 

r  •  l 


Cu.  — _ 

J  ar+  «'  «r+  i  (u>+0r) 

{' V,-iV, ) 


<FO)> 


v'12) 


«.+  »  (•«  +  Pr) 

From  this,  it  is  not  difficult  to  extend  the  above  identification  tech¬ 
nique  to  the  condition  under  which  only  output  (response)  data  exists. 


III.  IMPULSE  RESPONSE  METHOD 


The  Impulse  response  function  of  a  discrete  linear  steady  system 
with  N  degree  of  freedom  can  be  expressed  b y  the  inverse  Laplace  trans¬ 
formation  of  the  relation  between  the  transfer  function  and  the  com¬ 
plex  frequency  sr  and  the  complex  number  Ar  equation  (8):  36 


a„(  t )«  S 


(13) 
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p  and  q  are  the  excitation  and  measurement  coordinates  respectively. 


At  time  t  ,  let  and  note  .  The  impulse  res¬ 

ponse  at  time  t  can  be  expressed  as  (omitting  the  symbols  p,  q) 

where  n  is  the  power.  For  each  sampling  time  of  equal  difference  • 
i.(n«o,  1,  2,  ,  iV),  ,  we  can  write  (2  N  +  1)  eauations.  Let 

V  (1  c. ) 

n  ix-x,Hx-xn-  s  a*x*“  ° 

r  «1  "  “  0 

where  A'-#'4',  and  the  coefficients  ^(n = o ,  1,  2,  . .  2AO  ,  are  the 

autoregressive  coefficients  of  the  equation.  Once  these  coefficients 
are  obtained,  it  is  possible  to  obtain  the  complex  frequency 
*,=  —  a,+»0,  by  solving  the  complex  root  A-,— using  equation  (15) 

From  them,  we  can  obtain  the  modal  frequency  and  damp  ratio 
2,  . ,  N) 


3f,  . 


The  method  to  obtain  the  auto-regressive  coefficients  is  as 
follows:  multiply  both  sides  of  equation  (14)  by  an  and  let 


0  ,  1  ,  2  ,  . ,  2  N, 


After  rearrangement,  we  get 


(m  =*  0  ,  1,  2,  . ,  2N-  1)  Ut) 

n  *  0 

There  are  only  2N  independent  amplitude  values  in  the  impulse 
response  data.  Let  us  choose  a2NBl  to  standardize  it.  Then,  we  have 

iff 

2  “  A  .  (17) 


By  considering  random  interference  and  testing  errors,  the  total 
square  deviation  of  impulse  response  corresponding  to  M  sampling 
times  is 


m  /  tN  -  i  y 

/,( « )  -  2  2  a  (/-♦*) 

m  -  1  \  •  -  »  / 


The  extreme  value  {o)-(a„  «n  ®»*  . »  a*"-^T  corresponding  to  J^(a) 

min  can  be  obtained  from  the  following  set  of  equations 
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where 


2  R,.^.=  -R,.i»  («  =  o,  l,  2, . ,2N-l)  (15) 


/?,..=  S  hO—AhU  »♦•) 


is  the  correlation  function  of  impulse  response.  The  latter  is  an 
even  function  which  only  deviates  from  the  time  difference.  Equa¬ 
tion  (19)  can  be  written  in  the  following  matrix  form 


r  /?( o ) 

i?(l) . 

■-RC2N-  1  ) 

,  (  2iV )  s 

Af(l) 

/?(o) . 

■R(2N-  2  ) 

! 

;  -oi 

i  • 

/— s 

1 

CM  M1 

05 

1 

II 

R(2N  —  1  ) 

R  (2N  —  2)" 

0) 

1 

1  : 

(/?(!)  J 

The  least  square  method  is  as  follows:  Express  the  Impulse  res¬ 
ponse  in  the  form  of  sine  and  cosine  function  form  (notice  that 

~ar+  ,prt  A,  —  u,  +  ior) 

N 

A(/„)  =  2  ^  «'*''“(»rC08(Pr/«)-t;rsin(pr/_))  (21) 

r  •  1 

The  actual  measured  Impulse  response  value  at  tire  tm  is  denoted 
as  gm;  then  the  total  square  deviation  of  M  sampling  values  is 


where 


M  /  N  \» 

A(  “  *  *0  =  2  2  (c,mUr  -S,mVr)  ]  (22 

m  *  IV  r  -  x  J 


c,m= ■  2e”''-cos(Pr<«)»  2e"''«sin(pr/») 


The  extreme  value  ({/,  F)T“( «„  u, •••«„,  w,—vw)t  of  J2(u,v)  =  min 

can  be  obtained  from  the  following  set  of  linear  functions 


\a  "T w  l-r  ^  1 

U>T  bJL  V  J  L  Y  J 


where  A,  B,  and  D  are  Nth  order  matrices.  The  ith  row  jth  column 
elements  are  m  m 

All “  2  Crmfilmt  Bl{^  J]  *(«.*/■ 

at  ■  l  m  *  x 

M 

Du ™  2 
m  •  1 


,W.  •-  V 


X  and  Y  are  N  dimensional  vectors  whose  element  is  related  to  the 
measured  data 

M  M 

CimQwt I  K,  =  -  ^  '  SimQm 

in  ■  1  m  =  1 

The  original  data  of  impulse  response  system  identification  is 
the  impulse  response  tested  data  which  can  be  obtained  using  the 
following  method 

(1)  if  input  and  output  data  can  be  simultaneously  measured,  then 
the  frequency  response  can  be  obtained  first.  Then,  Fourier  inverse 
transform  can  be  used  to  obtain  the  impulse  response  data.  Both  can 
utilize  FFT  and  BSFA  techniques; 

(2)  if  the  input  (excitation  force)  cannot  be  directly  measured, 
then  a  random  decay  technique  can  be  used  which  triggers  the  record¬ 
ing  at  a  certain  voltage  with  respect  to  the  wide  frequency  band 
random  response  signal  of  the  system.  Then  it  is  followed  by  the 
sampled  data  total  average.  Assuming  that  the  forced  response  of 
the  voltage  triggering  sampling  signal  after  total  averaging  and  the 
response  caused  by  the  initial  triggering  condition  are  close  to  zero, 
the  obtained  characteristic  signal  can  be  described  by  the  impulse 
response . 

IV.  AUTOREGRESSIVE  MOVING  AVERAGE  MODEL  METHOD 

The  usual  expressions  of  a  dynamic  system  in  a  time  domain  are 
modal  functions  and  output  functions.  For  a  time  discretized  system, 
the  input/output  relation  can  be  expressed  by  the  difference  equation 

2N  J.V 

2  <Jix(  t  -  I  )  =  £  b.  v(  t  -  i  )  (  2^  ) 

«'  ■  0  I  «  0 

Assume  that  y(t)  is  an  independent  random  serial  input  with  an 
average  zero,  square  deviation  a2.  After  taking  random  measurement 
noise  into  account  without  losing  generality,  the  above  equation  can 
be  written  as: 


63 


'v  y- 


38 


w 


U*i!  -WUUJw*' 


iN  W 

y!  a,.xf  t  —  i  )—  b,y(  t  —  i) 

*— *  .  .  A 


(a„  =  60  =  1  ) 


(25) 


The  measured  response  value  of  a  certain  coordinate  of  the 
system  at  time  t  is 

f  (  t  )  =  —  a,x(  t  —  1  )  —  J;-*(  /  —  2  )  — —  /  —  2AT)  +>'(<)  (26) 

+  6,y(  t  —  1  )+6,y(  t  —  2  )  + -rk._sy(  t  —2 N) 

The  first  2N  terms  are  the  sum  of  the  measured  values  (auto- 
regressing)  before  time  t,  and  the  last  (2  N  +  1)  terms  are  the  mov¬ 
ing  average  of  the  input.  Therefore,  equation  (25)  is  also  called 
the  autoregressive  moving  (ARMA)  model.  The  autoregressive  coeffi¬ 
cients  a^  and  the  characteristic  values  *=—  have  the  following 

relation 

V  «»*'-  fl  CX-X.KX-XV  (27) 

i  «  0  i  *  1 

where 

After  identifying  the  autoregressive  coefficients,  it  is  possible  to 
solve  the  attenuation  coefficient  and  damping  free  vibrational 

frequency  (from  here  we  can  obtain  the  damp  ratio  and  modal  fre¬ 
quency).  The  moving  average  coefficients  bi,  however,  are  related 
to  the  vibration  model  information. 


The  coefficients  of  the  ARMA  model  can  be  identified  using  the 
maximum  approximation  method.  The  standard  function  can  be  chosen 
as  the  total  square  deviation  (approximation  function)  of  the  measure¬ 
ment  noise  M 

no,  6)«~ jf ^2  y(ty=cyy(0)  (2S) 

where  y(t)  is  the  measurement  noise,  cyy(O)  is  the  auto-correlation 
function  with  zero  time  difference,  and  M  is  the  sampling  data  points. 
It  is  possible  to  use  the  Newton-Raphson  method  to  solve  the  coeffi¬ 
cients  a b^  (i  =  1,2, — 2N)  to  be  identified  to  make  J(a,b)  =  min. 

If  the  input  cannot  be  measured,  cyy(O)  can  be  calculated  from  the 
auto-correlation  function  cxx(i)  of  the  response. 


The  advantages  of  maximum  approximation  parameter  estimation  are 
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approaching  without  deviation  (the  expected  estimation  value  is 
equal  to  the  real  value),  identical  intensity  (when  the  sampling 
point  increases,  the  estimation  value  converges  to  the  actual  value 
with  probability  equal  to  1),  and  statistically  effective  (the  least 
square  deviation  between  the  estimation  and  real  values  can  reach 
the  Cramer-Lao  lower  limit).  The  shortcoming  is  that  the  computation 
is  too  complicated  which  is  suitable  for  engineering  applications. 

Let  us  write  equation  (26)  (corresponding  to  M  sampling  time 
inoervals)  in  the  matrix  form 

<  JC  >  —  C/>3T{  9  >  -#-C  e  3  (29) 

where  (x)  is  the  vector  formed  by  the  sampled  response  values  at  t  = 

2N  +  1,  2N  +  2,  ....,  2N  +  M,  { e }  is  the  corresponding  measurement 

T 

error  vector,  (9)  =  (b^,  a2]^  the  vectors  to  be  iden- 

tified.  The  coefficient  matrix  is 

y(2N ),  y(2N+i) y(.2N+M-i ) 

—  x  (2jV),  —  x  (2iV+  1  ) . —  x  (2N  +  M—  1  ) 

y(D  y( 2) .  y(M) 

-*(1)  ~*(2) . -*(M) 

From  equation  (29),  the  least  square  estimation  of  the  parameter 
to  be  identified  can  be  obtained  from  the  following  equation: 

{0>=*CCP)CPDT)'‘CPDT<^>  (30) 

The  standard  square  deviation  between  the  estimated  value  and  the 
real  value  can  also  be  estimated 

CoKe,  0  )«£<(§-  0)(e- 0)T>— ;vyl2\v  < CCPHP^r'  (31) 

where  [P]  [P]  can  be  calculated  from  the  correlation  function  of 
the  testing  data. 

F (  I  ) . R( 2<V-  i  ) 

F(  0  ) . R  ( 2.V  —  2  ) 

1  )  R(2N-  2  ) . F(  0  ) 

cyy(  i  )  -cyx(  i  ) 

exy(  i  )  cxx(  i  )  „ 
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P(o) 

cpkp)t*  pc  i) 

.R(2N- 


V 


Under  the  condition  that  the  input  signal  (such  as  atmospheric 

turbulence,  gust)  cannot  be  observed,  the  least  square  method  can  be 

used  which  involves  the  structuring  of  an  autoregressive  (AR)  model 

first  and  then  deriving  the  auto-correlation  and  mutual  correlation 

functions  of  the  input  based  on  the  auto-correlation  function  of  the 

output  data.  Let  the  AR  model  be 

L 

S  ay*(t  -  (a,-l)  (32) 

where  y(x)  satisfies  '  ’ 8 

E{y(  O>-0i  E<*(  * 

ai  are  the  autoregressive  coefficients  to  be  identified  and  L  is  the 
order  to  be  determined.  Because  y  and  x  are  not  correlated,  we  find 
from  equation  (32)  that 

L 

2]  <*icxx(  k  —  i  )—  0  (  A  -  1 ,  2 ,  . ,  I) 


i  •  l 


(33) 

From  this  we  can  find  a^.  The  order  L  can  be  identified  from  the 
Akaike  information  guideline.  The  square  deviation  estimation  value 
of  measurement  error  is 


$*-  2  a<cx*(  • ) 
i  •  I 


(3M 


From  the  relation  between  input  and  output 

L 

*<«)-  2  f>(i)y(t-i ) 


(35) 


h(i)  is  the  impulse  response  function.  The  above  expression  is  also 
called  the  moving  average  (MA)  model.  From  equation  (35)  we  can 
derive  the  mutual  correlation  function 


exy( 


olh(  i  ) 
0 


»  >  0 

i  <0 


(36) 


h(i)  can  be  obtained  from  the  following  recurrence  formula 


fc(0)-l,  /»(<)-“  £  a(i)fc(<-»)  (1-1,  2,  •••) 


(37) 


»  *  1 


From  this,  the  original  data  needed  to  identify  the  parameters 
and  the  square  deviations  of  the  real  values  of  the  ARMA  model  by 
equations  (30)  and  (31)  can  be  obtained  from  the  sampled  data. 
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V.  CONCLUSIONS 


1.  Based  on  the  frequency  response  data  obtained  in  ground 
vibration  measurements,  it  is  possible  to  obtain  sufficiently  accur¬ 
ate  structural  modal  coefficients  using  the  transfer  function  identi¬ 
fication  method.  It  has  been  a  success  for  the  flutter  model  of  a 
wing  with  external  bodies.  This  technique  can  be  extended  to  the 
testing  of  the  entire  aircraft  or  to  determine  the  pure  modal  excita¬ 
tion  shaker  forces  of  multiple  points  through  a  transfer  function. 

In  the  wind  tunnel  and  flight  flutter  tests,  it  is  possible  to 
obtain  the  original  data  using  random  or  transient  excitation  PFT 
analysis,  including  the  situation  with  response  excitation  alone. 

In  order  to  verify  and  modify  the  finite  element  structure  dyna¬ 
mic  modal  mathematical  model,  more  reliable  and  accurate  complete 
modal  parameters  are  required.  We  suggest  using  steady  state  sinus¬ 
oidal  excitation,  tracking  filtering  to  analyze  the  frequency  response 
data  for  identification. 

The  transfer  function  method  can  also  obtain  the  distribution 
of  zero  and  extremum  points  based  on  the  aerodynamic  elasticity  cal¬ 
culation  of  the  aircraft  or  the  rational  transfer  function  model 
formed  according  to  the  experimental  data  in  order  to  facilitate  the 
analysis  and  synthesis  of  the  flutter  of  the  stabilizer  system  and 
the  flutter  active  control  system. 

2.  The  impulse  response  method,  in  addition  to  the  fact  that  it 
is  possible  to  use  the  inverse  Fourier  transform  of  the  measured 
data  as  the  original  data,  is  especially  suitable  for  the  condition 
with  response  data  only.  At  this  time,  the  random  attenuation  method 
can  be  used  to  obtain  the  original  data  without  any  data  analyzer. 

It  is  suited  for  the  Identification  in  wind  tunnel  and  flight  tests. 

3.  The  autoregressive  moving  average  model  method  is  an  effect¬ 
ive  statistic  parameter  estimation  technique.  Because  the  random 
noise  effect  Is  confronted  face-to-face  in  this  method,  and  also 
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because  a  double  smoothing  correlation — least  square  method  is  used 
it  is  expected  to  obtain  higher  identification  accuracy  with  strong 
interference.  This  method  can  also  simultaneously  obtain  the  stand¬ 
ard  square  deviation  between  the  real  and  the  estimated  values  to 
further  ensure  the  reliability  of  the  results.  Furthermore,  this 
method  can  be  directly  based  on  sampled  data  without  the  use  of  spe¬ 
cial  analytical  instruments. 
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SYSTEM  IDENTIFICATION  AND  AIRCRAFT  FLUTTER 


Zhang  Lingmi 

(.Nanjing  Aeronautical  Institute ) 

Abstract 

After  reviewing  the  application  of  system  identification  to  aircraft  flutter 
research«  three  methods  of  system  identification  are  presented  in  this  paper. 
Emphases  are  put  on  the  main  points  of  identification  algorithm  and  its  theo¬ 
retical  basis.  The  problems  related  to  data  processing  are  also  discussed.  The 
Transfer  Function  Method  based  on  complex  modal  analysis  and  optimization 
technique  can  identify  all  of  the  modal  parameters  accurately.  Furthermore, 
this  method  has  been  applied  to  the  flutter  model  test  of  a  wing  with  external 
bodies  successfully.  The  method  can  be  extended  to  the  cases  of  wind  tunnel 
and  flight  tests  which  provide  response  data  only.  Adopting  the  correlation 
and  least  square  technique.  Impulse  Function  Method  and  Autoregressive  Mo¬ 
ving  Average  Method  can  gain  considerable  accuracy  in  identification  of  measu¬ 
rement  data  conteminated  with  rather  strong  noise  interference.  In  no  need 
for  a  special  signal  analysis  instrument,  it  is  possible  to  make  direct  use  of 
measurement  sampled  data  from  artificial  random  or  transient  excitation  and 
response  data  from  natural  random  excitation  with  these  two  time  domain 
methods. 
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APPLICATION  OF  MULTIPLE  DYNAMIC  ABSORBERS  TO  REDUCE 
THE  VIBRATION  LEVEL  OF  A  COMPLEX  CANTILEVER  STRUCTURE 


TIan  Qlanli 

Institute  of  Mechanics,  Chinese  Academy  of  Sciences 


ABSTRACT 

A  complex  cantilever  structure  has  two  closely  positioned 
resonance  peaks  near  20  Hz  causing  the  destruction  of  the 
root  of  the  structure.  This  paper  presents  a  method  to 
use  multiple  dynamic  absorbers  to  reduce  its  vibration 
level.  In  order  to  overcome  the  drawback  of  the  usual 
tuned  absorbers  which  are  very  sensitive  to  the  structural 
frequency,  six  absorbers  were  hung  in  a  given  section  of 
the  structure.  The  absorbers  can  vibrate  along  any  direct¬ 
ion  with  the  structure.  Furthermore,  the  stiffness  and 
damping  parameters  of  these  six  absorbers  could  be  different 
from  one  another  so  that  the  absorbing  frequency  range  could 
be  widened.  In  order  to  find  the  optimum  parameters  and  to 
study  the  effect  of  parameter  variation  on  the  structure 
response,  large  amounts  of  response  curves  were  calculated. 
Because  the  structural  damping  is  very  small  and  the  absor¬ 
ber  damping  is  large,  therefore,  it  is  a  non-proportional 
damping  dynamic  analysis  problem.  In  this  paper,  this  prob¬ 
lem  was  solved  by  an  eigen  solution  method  and  a  modal  syn¬ 
thesis  method. 


I .  INTRODUCTION 


For  a  cantilever  structure  in  resonance  at  the  base  frequency, 
a  large  stress  is  produced  at  its  root  which  easily  leads  to  destruc¬ 
tion.  Use  of  dynamic  absorbers  can  solve  this  problem  very  well. 
However,  for  a  complex  cantilever  structure,  there  are  usually  several 
resonance  peaks  in  the  vicinity  of  the  base  frequency.  In  addition, 
due  to  the  fact  that  the  direction  of  excitation  vibration  is  not 
specified,  these  resonance  peaks  are  frequently  coupled.  The  fre¬ 
quency  parameters  in  the  long  period  of  the  working  process  of  the 
cantilever  and  absorbers  may  deviate  from  the  original  designed  numer¬ 
ical  values.  Hence,  it  is  impossible  to  use  classical  tuned  absorbers 
to  treat  this  problem.  We  adopted  the  method  of  installing  multiple 
absorbers  with  different  parameters  on  a  section  of  the  structure  in 
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order  to  widen  its  absorbing  frequency  band  and  to  reduce  the  vibra¬ 
tion  level  in  any  direction  by  vibrating  in  any  direction  following 
the  structure. 


We  all  know  very  well  that  the  commonly  used  finite  element 
analysis  programs  such  as  SAP  can  only  calculate  the  dynamic  prob¬ 
lems  with  proportional  damping.  Here,  our  damping  is  concentrated 
on  the  absorbers  and  these  absorbers  are  hung  on  a  section  of  the 
structure  which  is  apparently  a  non-proportional  damping  problem. 
Therefore,  we  used  a  complex  eigen  solution  method  and  a  modal  syn¬ 
thesis  method. 

II.  USING  A  MODAL  SYNTHESIS  METHOD  TO  CALCULATE  THE 

DYNAMIC  RESPONSE 

The  special  feature  of  the  modal  synthesis  method  is  that  the 
modal  solution  of  each  substructure  is  obtained  individually  first. 
Then  the  eigen  vectors  of  the  substructures  are  synthesized  into  the 
Ritz  vector  of  the  structure  in  order  to  solve  the  eigen  solution  of 
the  system  in  the  subspace.  Here,  the  damping  of  the  structure 
itself  is  very  small,  while  the  damping  of  the  absorbers  is  larger. 
Therefore,  we  divided  them  into  two  substructures.  The  dampless 
mode  of  the  structure  itself  is  solved  by  the  SAPN  program.  A  few 
lower  order  modes  after  the  cutoff  are  combined  with  the  concentrat¬ 
ing  parameters  of  the  absorber  to  carry  out  the  synthesis.  Now,  the 
first  few  terms  of  the  natural  frequencies  calculated  by  the  SAPN 
program  and  the  experimental  results  are  listed  in  Table  1  for 
comparison . 


Table  1  The  few  terms  of  natural  frequencies  of  a  cantilever  structure 


frequency 

order 

/. 

/» 

ft 

/« 

/. 

ft 

calculated  value 

20.29 

20. S3 

SI.  3 

14S.0 

141. S 

258.3 

sinusoidal. test 

20.19 

21.42 

IS 

147 

148 

240 

random  test 

20.22 

21.41 

13.10 

190.41 

191.23  ! 

299. S7 

Because  f^  and  f2  are  very  closely  connected,  they  correspond 
to  the  bending  vibrations  in  two  main  directions,  respectively. 


A 


In  order  to  demodulate  them,  we  used  the  mechanical  resistance  test¬ 
ing  method.  In  the  20-22  Hz  frequency  band,  we  used  slow  speed 
scanning  and  the  phase  and  amplitude  were  recorded  at  high  speed  us¬ 
ing  an  X-Y  recorder.  Then  the  resistance  circle  was  plotted  manually 
point-by-point.  Therefore,  in  this  frequency  band,  the  frequency 
resolution  could  reach  0.01  Hz.  The  random  excitation  vibration  was 
controlled  by  B  *K  3380.  The  response  was  recorded  on  magnetic  tape 
and  spectral  analysis  was  carried  out  on  C-F  700. 

The  actual  measured  base  frequency  vibration  model  basically 
agrees  with  the  calculated  results.  Therefore,  we  began  the  synthesis 
using  the  calculated  modal  parameters  and  the  absorbers.  As  shown 
in  Figure  1,  after  installing  a  dynamic  absorber  at  point  z  of  the 
cantilever  structure,  it  is  equivalent  to  the  exertion  of  a  concen¬ 
trated  force  at  that  point 

(1) 

Here  */.« ■>*/.»(  1  +*P)  is  the  complex  stiffness  of  the  complex  spring  of 
the  absorber,  8  is  the  damping  coefficient,  i  —V-\ . 

In  the  finite  element  calculation,  it  is  necessary  to  discretize 
the  structure  into  several  nodal  points  and  each  nodal  point  can  have 
0-6  degrees  of  freedom.  If  the  beam  element  is  one  which  neglects 
shear  and  rotational  inertia,  then  each  point  has  two  degrees  of 
freedom.  Therefore,  in  the  following  we  will  use  the  jth  degree  of 
freedom  to  represent  the  Zj  point  mentioned  above.  By  assuming  the 
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total  degrees  of  freedom  after  structural  discretization  to  be  N, 
then  the  displacement  vector  X  is  an  N  dimensional  vector.  In  a 
structure  with  m  degrees  of  freedom  with  p  absorbers  when 

the  structure  is  under  motion  with  a  base  acceleration  equal  to  x, 
the  equation  of  motion  of  the  system  is 

(2) 

M*+KX=M(/I*i+F 

where 


F~ 


*:.I  (y:.t 


-x.) 


I 


(3) 


The  subscript  j  represents  the  jth  degree  of  freedom,  the  subscript 
1  represents  the  1th  absorber,  the  superscript  *  represents  the  com¬ 
plex  number  of  the  quantity,  { 1)T-  (  l ,  l  •••  i )  .  Because  we  installed 
six  absorbers  of  various  parameters  on  a  section  of  a  structure, 
therefore,  in  equation  (3)  corresponding  to  the  same  x^  ,  it  is  poss¬ 
ible  to  have  1  number  of  y^  In  order  to  unify  the  two  in  one  coor¬ 

dinate  system,  let  us  introduce  a  matrix  3 


1  0  0 
0  1  0 
0  1  0 

:  0  1 
:  0  1 

!  :  0 


(4) 


8  is  a  p  x  m  matrix  in  which  the  number  of  rows  corresponds  to  the 
number  of  absorbers  and  the  number  of  columns  is  equal  to  the  degrees 
of  freedom  of  the  added  absorbers.  Its  element  corresponding  to  x^ 
in  each  column  is  equal  to  1,  and  the  remaining  elements  are  zero. 

Let 


(5) 


T  is  a  p  x  p  diagonal  matrix  and  its  elements  are  the  complex  stiff¬ 
ness  of  the  absorbers.  Let  y  be  a  p  dimensional  vector  whose  elements 

are  the  displacements  y,  .  of  the  absorbers.  The  sequence  corresponds 

J  » 1  * 

to  the  diagonal  elements  of  T  .  Therefore,  the  F  in  equation  (3) 
can  be  written  as 


F“PTT#  {y—fiX/} 


(6) 


„  A 


HU 
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Xj  is  the  displacement  vector  corresponding  to  the  part  of  degrees 
of  freedom  of  the  attached  absorber  in  the  total  displacement  vec¬ 
tor  of  the  structure  X.  It  is  related  to  X  through  the  projection 
matrix  P: 

X,-PX  (7) 

r  o  i  o  o . 1 


oooi . 

oooooi 


P  is  an  m  x  N  matrix  [1].  Only  the  j^  element  in  each  row  which 
corresponds  to  the  degree  of  freedom  of  the  attached  absorber  is  equal 
to  1,  and  other  elements  are  zero. 


Prom  equation  (6),  we  can  see  that  P  is  an  m  dimensional  vector. 
T 

Multiply  it  by  P  from  the  front,  so  that  it  is  extended  into  an  N 
dimensional  vector.  Hence,  equation  (2)  can  be  written  as 

Mx+KX+PJfiJTmtPX-PTtJTmy- ( 9 ) 

Simultaneously,  the  equation  of  motion  p  absorbers  can  be  written  as 

jir +t*y-t*ppx- -am*.  (10) 


where 


Is  a  p  x  p  diagonal  matrix 


Ri,i  is  the  mass  of  the  1  absorber  attached  to  the  (JU"  degree  of 
freedom  of  the  structure.  Combining  equations  (9)  and  (10),  we  get 
the  equation  of  motion  of  the  synthesized  structure: 

I'M:  O')/*)  ,  fK  +  PTPJT*pPi  -PtPtT*  1  fX ) [M{\  )).. 

[  0*7  R J  liri  + 1 — T*jJP | . T* j  lY  J  ( R  <  1  )  i  •  (12) 

In  equation  (12),  T  is  a  complex  number  diagonal  matrix.  Therefore, 
it  is  a  (N  +  p)th  order  complex  number  matrix  equation.  The  order  N 
of  a  complicated  structure  is  very  high,  and  it  takes  a  great  deal  of 
computer  time  to  solve  the  complex  eigen  solution.  For  this,  let  us 
take  the  first  n  terms  of  modes  of  the  structure  as  the  Hitz  vector 
of  the  synthesis  system.  Let  x  =  Vq  and  substitute  it  only  into 
equation  (9)  where  V  is  the  matrix  formed  by  the  first  n  terms  of  the 
eigen  vectors  of  the  structure.  Again,  multiply  it  by  VT  from  the 
front.  We  notice  that 


VTMV=/  (I  is  the  unity  matrix) 


vtkv-a 


t  h 

is  the  r  order  eigen  value 
of  the  structure 


Then,  we  have 

/«+Cy\  +  VTPTpTT*PPVDq-VTPTpTr*Y--VTM{  1  }3e,  (13) 

Let  PV-V,  PTT*-T!  PtT*P-TJ 

Then  equation  (12)  becomes 


a+v,tt?v, 
. -tTv," 


i-V/TT 

:  T* 


]{$}--{ 


VTAf {  1 } 

p  { 1  > . r* 


(14 ) 


Through  the  aforementioned  cutoff  modal  transformation,  the 
order  of  the  square  matrix  decreases  from  (N  +  p)  to  (n  +  p)  and 
usually  n<<N.  Therefore,  the  computational  load  is  greatly  reduced. 
This  point  is  especially  significant  in  dealing  with  the  absorber 
problem.  Because  the  base  frequency  component  occupies  an  important 
portion  in  the  cantilever  structure  vibration  response,  it  is  only 
necessary  to  keep  the  first  few  orders  of  modes  to  obtain  a  solution 
of  sufficient  accuracy.  Here,  our  structure  has  54  degrees  of  free¬ 
dom  after  discretization,  i.e.,  N  =  5**.  The  SAP  program  is  used  to 
obtain  the  first  10  orders  of  eigen  values  using  an  iteration  method 
in  the  substructure  space.  The  sixth,  third  and  second  order  eigen 
vectors  were  selected  to  form  the  vector  V  to  be  substituted  into  the 
left  side  of  equation  (14).  Use  Q.  R.  Yugenmoyacoby : s  complex  eigen 
solution  program  [2]  to  calculate  the  complex  eigen  values  and  complex 
eigen  vectors,  and  the  results  obtained  indicate  that  they  are  almost 
identical  to  those  obtained  when  n  =  6,  3,  2.  The  difference  is 
smaller  than  one  one-thousandth. 


III.  COMPLEX  RESPONSE  FUNCTION 


Equation  (14)  is  a  complex  number  matrix  equation;  let  Q 
From  the  complex  eigen  solution  program,  we  can  obtain  the  eigen 
solution  of  the  homogeneous  equations  on  the  left  end  of  equation  (14). 


Q-$*Z 


(15) 
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•*  is  the  complex  eigen  vector  matrix. 


Substituting  equation  (15)  into  equation  (14)  and  multiplying 
it  by  0*T  from  the  front,  we  can  obtain  the  following  using  ortho¬ 
gonality  and  normalization  characteristics: 

z,(  t  )+°>K  i  +»n,)2,(  /  )  =  -#?T[  j1  '  }*•  (16) 

r  -  1  .  2  , . « 

$*T  is  the  r^  row  in  ®*T 


Let 


*r{  ,-d* 

(17) 

Conducting  Fourier 

transform  on  equation  (16),  we  get 

Here 

Zr(a>)  =  H,(<o)P,(a>) 

(18) 

(19) 

(20) 

(21) 

The  denominator  of  Hr(“>)has  two  roots,  ± (■»,  +  (  1  +  s'n,)''‘ 

where  -<n,(  l  +in,),/‘  is  located  below  the  real  number  axis  of  the  com¬ 
plex  plane.  Therefore,  «>(“)  is  unstable  in  the  frequency  region. 

For  this  purpose,  we  adopted  the  suggestion  of  [3]  to  use  the  sum  of 
a  pair  of  complex  conjugate  values  to  represent  the  real  solution, 
i . e . ,  to  let 

Zr(  0>  )  -  H,(  <0  )p,(  CO  )!.>,  +  Hr(  <0  )  >,(  ®  )Lo  (22) 

Here  5,(<o),  p,(a>)  were  taken  from  equations  (21)  and  (22),  respectively. 
P*(®)  are  their  complex  conjugates. 

ft 

Let  us  divide  the  $  matrix  in  equation  (15)  into  two  parts  which 
correspond  to  f-9l  ,  respectively,  i.e.. 


r in*-;?! 

then  {q)*=Q\Z  (23) 

Hence,  the  displacement  vector  of  the  structure  is 

x(ci>)=vc<J>;z(co)  L>,+iiz(“)  l.<0 

IV.  RANDOM  VIBRATION  RESPONSE 


The  Fourier  transform  of  the  bending  moment  of  the  root  of  a 
cantilever  structure  during  random  vibration  is 


MO)-  {hVMAVX O) 


Here  (MT  is  the  transport  of  the  vector  formed  by  the  root 
heights  of  all  the  nodal  points  after  the  discretization  of  the  struc- 

i.  y. 

ture.  V  is  the  rc  column  vector  in  the  V  matrix,  and  <J>*,  is  the 

y  u  *  t  h 

r  row  nr  column  element  in  the  matrix 


Let 

o!(h)TMV,=»c, 

(26) 

then 

Bmf  Bm 

r  r 

(27) 

M(  CO  )  -  2  (BlZmi  <0  )|#>f  +  BlZmi  <*>  )!*<,) 

tarn 

(28) 

Its  complex  conjugate 

•m 

is 

M  ( » )  -  2(  Btz-(  ®  )!.>,+ Biz. (  ®  )!.<•) 

(29) 

The  spectral  density  of  the  bending  moment  at  the  root  of  the 
structure  is  S„ 0)=^Lim  *  S  2 CbSbI5J62h.hJ«>« 

+  BlBiblDlH.H «> )  ^  30 ) 

Let  +  (3D 

then  BlBlblbi-Am,.-iBm..  (32) 


Substituting  it  into  equation  (30)  and  integrating  with  respect 
to  a)  from  -x  to  +«,  the  mean  square  value  of  the  root  bending  moment 
can  be  obtained  as 
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If  iS;-,v,(®)  varies  slowly  in  the  frequency  band  of  interest, 
be  approximated  by  the  white  noise  spectrum  S.  Then 

(  I  2a>.(ai  +  6i)t<»i(  1  +>li)+^(  1  +  hi) -  2o>i.®i(  1  -h-h.))  J 


it  may 

(3*0 


where 


a. 


=  i  VT+TI+  i 

v'  2  \ 


When  ®«,  ®.  are  separated  relatively  far  apart,  and  TU~'n.<0.1 
the  coupled  term  of  m  and  n  can  be  neglected.  This  caused  the  approx¬ 
imate  expression  in  equation  (2.130)  of  [4]. 


V.  APPLICATIONS 


From  equations  (21)  and  (35),  we  can  see  that  the  response  func¬ 
tion  and  the  mean  square  value  of  response  are  inversely  proportional 
to  the  modal  damping  regardless,  whether  it  is  simple  harmonic  vibra¬ 
tion  or  random  vibration.  Of  course,  it  is  also  related  to  the  vibra¬ 
tional  model  parameters.  Basically,  the  stress  is  proportional  to 
the  square  of  vibration  model.  By  adjusting  the  weight,  position  and 
parameters  of  the  absorbers,  It  Is  possible  to  change  the  eigen  solu¬ 
tion  of  the  system;  i.e.,  it  might  change  the  response  value.  Here, 
due  to  the  limitation  in  structure,  the  weights  and  positions  of  ab¬ 
sorbers  are  defined.  We  are  only  able  to  change  their  stiffness  and 
damping  parameters.  Because  the  base  frequency  of  a  cantilever  struc¬ 
ture  has  two  major  directions  of  bending  vibration,  therefore,  the 
absorbers  are  divided  into  two  groups.  Each  group  has  three  absorbers 
corresponding  to  one  major  bending  direction.  Their  stiffness  and 
damping  may  be  different.  The  stiffness  coefficient  k,  ,  ranges  from 
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0. 8-3.0  kg/cm  and  the  damping  coefficient  B  varies  in  the  range  of 
0.05-1.0.  With  regard  to  various  complex  stiffness  combinations,  a 
large  amount  of  numerical  calculations  has  been  carried  out.  We  dis¬ 
covered  that  when  the  stiffness  coefficients  of  the  three  absorbers 
are  1.5,  2.0  and  2.5  kg/cm,  respectively,  and  when  the  damping  coeffi¬ 
cients  are  0.2,  the  vibration  absorbing  effectiveness  is  optimal. 

Now,  a  comparison  of  these  optimal  response  curves  and  another  series 
of  response  curves  obtained  with  identical  stiffness  coefficients 
and  S  equal  to  1.0  is  shown  in  Figure  2. 

It  is  worthwhile  mentioning  that  when  this  structure  was  calcul¬ 
ated  based  on  an  equivalent  single  degree  of  freedom  body  system,  the 
optimal  stiffness  coefficient  is  1.6^4-1.66  kg/cm  and  the  optimal 
damping  coefficient  is  0.19.  However,  due  to  the  addition  of  six 
absorbers  to  reduce  the  vibration  along  the  two  major  bending  direct¬ 
ions,  the  calculation  formulas  single  degree  of  freedom  system  are 
no  longer  applicable.  However,  calculated  results  showed  that  a  damp¬ 
ing  value  of  about  0.2  is  still  proper.  Figure  2  explains  that  the 
effectiveness  is  reduced  when  the  damping  value  is  too  large.  When 
the  frequency  of  the  absorber  is  lower  than  the  base  frequency  of 
the  structure,  the  absorber  consumes  energy  during  resonance  of  the 
structure.  [5]  had  used  this  property  to  solve  the  vibration  problem 
of  the  SMS  aircraft.  In  that  case,  the  damping  of  the  absorbers 
must  be  large.  However,  the  effect  was  not  as  good  as  the  multi¬ 
dimensional  absorbers  because  there  is  always  one  absorber  resonating 
when  the  structural  frequency  changes.  Both  calculated  and  experi¬ 
mental  results  showed  that  when  the  structural  frequency  and  absorber 
frequency  vary  by  305,  the  shock  absorbing  effect  can  still  reach  505. 

The  random  tests  and  sinusoidal  frequency  scanningtest  in  this 
paper  were  completed  by  comrades  Li  Yen  ping  and  Wang  Danfung.  The 
decomposition  of  the  base  frequency  was  performed  by  the  author 
himself.  The  calculations  were  carried  out  by  comrades  Li  Shen  zhang 
and  Liu  Dekong. 
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Fig.  2  The  frequeacy  response  curves  o{  a  structure  with  absorbers  having 
diUreat  parameters 
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APPLICATION  OF  MULTIPLE  DYNAMIC  ABSORBERS  TO 
REDUCING  THE  VIBRATION  LEVEL  OF  A  COMPLEX 
CANTILEVER  STRUCTURE 


Tian  Qianli 

(Institute  of  Mechanics ,  Chinese  Academy  of  Sciences') 

Abstract 

In  a  complex  cantilever  structure  there  are  some  closed  resonance  peaks  in 
the  vicinity  of  20Hz>  causing  a  serious  bending  moment  at  its  root.  The  appli¬ 
cation  of  multiple  dynamic  absorbers  to  reducing  its  vibration  level  is  propos¬ 
ed  in  this  paper.  Six  absorbers  are  hung  on  a  given  section  of  the  structure 
to  overcome  the  drawback  of  the  usual  tuned  absorbers,  i.  e.  excessive  sensi¬ 
tivity  to  the  tuning  parameters.  They  can  vibrate  in  all  directions  following 
the  structure,  but  their  stiffness  and  damping  parameters  of  these  absorbers 
are  different  from  each  other,  so  that  their  frequency  range  is  made  wide  enou¬ 
gh  to  cover  the  resonance  frequencies. 

In  search  of  the  conditions  for  minimizing  the  bending  stress  of  the  struc¬ 
ture  and  for  the  sake  of  studying  the  effects  of  the  parameters  on  the  dynamic 
response,  a  great  number  of  response  curves  at  the  top  of  the  structure,  bearing 
the  harmonic  excitation  from  the  base  movement,  are  calculated  as  the  pa¬ 
rameters  of  these  absorbers  vary  in  a  considerable  range.  Since  the  damping 
of  the  structure  is  very  small  and  that  of  absorbers  are  large  enough,  so  it  is 
a  dynamic  analysis  problem  with  non-proportional  damping.  This  problem  has 
been  solved  by  a  complex  eigen-solution  method  and  a  modal  synthesis  method 
in  the  present  paper. 
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A  SPECTRAL  APPROACH  FOR  ANALYZING  THE  VIBRATION  OF  A 
PERIODIC  STRUCTURE  WITH  RANDOM  PARAMETERS 


Huang  Wen  hu  * 

(Harbin  Institute  of  Technology) 


ABSTRACT 

In  order  to  explore  the  effect  of  fabrication 
deviation  of  the  blades  in  a  circumferentially  closed 
turbine  blade  assembly  on  the  vibration  of  the  turbine 
blades,  this  paper  used  a  structure  model  with  periodic 
random  parameters  as  an  approximation  of  the  blade  struc¬ 
ture.  In  addition,  a  spectral  method  was  presented  to  anal¬ 
yze  the  vibration  of  this  structure.  Assuming  the  stand¬ 
ard  deviations  of  the  structural  parameters  are  small, 
therefore,  it  is  possible  to  use  a  perturbation  method. 

The  periodic  random  structural  parameters  are  expanded 
into  Fourier  series,  so  that  the  free  vibration  and 
forced  vibration  of  the  structure  can  be  solved.  Then, 
the  frequency,  vibration  mode,  resonance  amplitude  and 
square  deviation  estimate  can  be  obtained.  The  orthogon¬ 
ality  of  the  main  vibration  modes  were  proven.  The  spec¬ 
ial  conditions  for  resonance  of  this  structure  were  anal¬ 
yzed.  The  examples  showed  that  the  analyzed  results  and 
experimental  results  have  the  same  order  of  magnitude. 


MAJOR  SYMBOLS 


At  Bt  a,  b,  e,  d - amplitude  coefficients 

/  (  O ,  /  ) - excited  vibration  force 

F(  e ) —  spatial  function  of  the  exciting  vibration  force 


h - number  of  harmonics  of  exciting  force 

l - order  of  the  Fourier  series 

m,  n ,  r - order  of  vibrational  mode  (nodal  diameter  no.  of  vibration  mode) 

tn.(  0  ),  k.(  0  ),  cg(  0  ),  d.(  0  )-  distributed  mass,  stiffness,  stiffness  of  connect- 

.  inq  and  damping 

P(  0  ),  Q  (  0  )— i local  frequencies  of  blade  and  connecting  part 

Pit  - average  val  ues  of  aforementioned  1  ocal  frequencies 

p(  q  )  Q(  Q  )  Af (  0  ) _ random  functions  of  distribution  stiffness,  connecting 

0  a  a  -standard  deviations  of  above  .  stiffness,  mass  parameters 

awi“UQlu  usviauwB  random  functions 

f  ,  g,  h — Courier  coefficients  of  above  random  functions 


x  (  0 ,  t  ) — displacement 
fc ,  ’l  >  t  —  smal  1  parameters 

Af(0)“i'(0)  +  4*(  0  )  +  nu(  0)4-  tui(  0  ) - vibration  mode  function 

X*+t»1*+Tlv*+tP* -  natural  frequency 


hQ — exciting  force  frequency 
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a,  0.  Y,  & - phase  angles 

e - damping  coefficient 

®  >  *  spatial  coordinates 
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I.  INTRODUCTION 

The  blades  of  a  turbine  are  usually  connected  into  a  group. 
Sometimes  through  the  use  of  various  types  of  connecting  parts,  such 
as  black  crows,  belts,  shoulders,  elastics,  it  is  possible  to  circum¬ 
ferentially  close  the  blades  of  the  turbine  disk  to  form  a  group. 

This  type  of  circumferentially  closely  connected  blade  group  has  cer¬ 
tain  advantages  In  avoiding  resonance.  With  regard  to  the  vibration 
of  such  a  structure,  references  [1-31  presented  a  method  to  calculate 
its  vibrational  characteristics  under  the  assumption  that  all  the 
blades  had  the  same  vibrational  characteristics.  In  addition,  a 
"triple  point"  condition  for  creating  resonance  on  the  rotating  disk 
turbine  blade  group  was  also  proposed.  This  means  that  resonance 
only  occurs  when  (1)  the  exciting  force  frequency  kfl  is  equal  to  the 
natural  frequency  wm  of  the  group  of  blades  (i.e.,  kft  =  u>m)  and  (2) 
the  excitating  force  harmonic  number  k  is  equal  to  the  nodal  diameter 
number  m  of  the  natural  vibration  mode  of  the  blade  group  (i.e.,  k  = 
m).  This  conclusion  could  only  be  obtained  under  the  assumption  that 
the  mechanical  properties  of  all  the  blades  and  connectors  were  identi¬ 
cal.  However,  in  reality,  It  Is  unavoidable  to  have  fabrication  devia¬ 
tions  between  blades  and  their  mechanical  properties  cannot  be  com¬ 
pletely  Identical.  Therefore,  the  above  conclusion  needs  some  correct¬ 
ion.  The  purpose  of  this  paper  Is  to  explore  the  effect  of  fabrica¬ 
tion  deviations  of  the  blades  on  the  vibration  of  the  blade  group. 

In  order  to  simplify  the  analysis,  we  chose  the  condition  under 
which  the  blade  distribution  on  the  turbine  was  very  dense  as  the 
limiting  case.  The  closed  loop  with  lateral  spring  support  was  select¬ 
ed  as  an  approximate  mechanical  model  to  obtain  a  differential  equa¬ 
tion  with  random  parameters. 
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In  the  random  vibration  problem,  when  the  parameters  on  the 
left  side  of  the  differential  equation  are  defined  quantities  and 
only  the  right  side  is  a  random  function  of  time,  this  type  of  prob¬ 
lem  has  been  studied  extensively  in  the  literature  with  matured 
results.  However,  with  regard  to  a  dynamic  structure  system  with 
random  parameters,  there  is  not  toe  much  available  in  the  litera¬ 
ture.  Furthermore,  there  is  a  lack  of  a  general  treatment  method. 
References  [M]  and  [5]  used  a  perturbation  method  to  study  the  free 
vibration  of  a  column  beam  with  random  structural  parameters.  The 
self-correlation  and  mutual  correlation  functions  of  random  structur¬ 
al  parameters  in  this  paper  were  given  manually  in  an  exponential 
function  form. 


This  paper  discusses  the  free  and  forced  vibrations  of  the  afore 
mentioned  closed  loop.  First,  by  assuming  that  the  standard  devia¬ 
tions  of  the  random  structural  parameters  are  infinitesimal  quanti¬ 
ties,  a  perturbation  method  is  applied  to  obtain  the  solutions.  By 
considering  that  the  closed  structure  treated  is  a  periodic  structure 
this  paper  proposes  to  use  a  spectral  method  to  find  the  solutions. 
The  random  structural  parameters  are  expanded  into  Fourier  series  to 
enable  the  solution  of  the  differential  equation  to  be  expressed  in 
terms  of  the  Fourier  coefficients  of  the  structural  parameters  in 
order  to  obtain  the  natural  frequencies,  natural  vibration  modes, 
resonance  amplitude  and  its  square  deviation  estimation  of  the  struc¬ 
ture.  The  results  of  examples  showed  that  the  analyzed  results  were 
in  the  same  order  of  magnitude  as  the  experimental  results  obtained 
in  [6,7]. 

II.  BASIC  EQUATIONS 


The  mechanical  nodel  discussed  In  this  paper  is  a  closed  loop 
with  lateral  spring  support  (see  Figure  1).  The  differential  equa¬ 
tion  of  motion  is 


+  k^e— 


d 

<)0 


~/(e,  t ) 


the  boundary  condition  is 


(1) 


62 


Fig.  1  Mechanical  model 


(_*,  ()-*(*,  i),  n  "’*0  l(*.  '> 


the  initial  condition  is 


Let  us  assume  that  the  .system  undergoes  a  harmonic  oscillation 
under  a  harmonic  exciting  force.  Let  mXQ  be  the  average  value  of 
the  mass  function;  and  the  exciting  force  be 

(3) 

/(e,  />-w.,F(ey* 

The  solution  to  the  differential  equation  (1)  is 

*(8,  o-xcey*  (3’) 

By  introducing  the  symbols: 


m«„  ”*««  '"•» 


then  we  can  obtain  the  differential  equation  of  the  vibration  mode 
from  the  partial  differential  equation  (1)  as 


The  boundary  condition  is 


III.  THE  RANDOM  STRUCTURAL  PARAMETERS 


Due  to  the  fabrication  deviation  which  exists  in  reality,  the 
structural  parameters  m  (0),  Ic  ( 9 )  and  c  (9)  are  generally  not  con- 
stants.  Instead,  they  are  random  functions  oscillating  around  some 
constant  values.  Also  because  what  is  discussed  here  is  a  closed 
structure,  these  functions  are  also  periodic  functions.  Parameters 
p*(6)  and  q2(0)  are  also  periodic  random  functions  oscillating  around 
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their  average  values  p2  and  q2  .  For  simplicity  in  this  paper  we 
assumed  that  the  damping  is  a  constant. 


Let  us  assume  that  the  scatter  of  a  structural  parameter  is  an 
infinitesimal  quantity.  Then  these  parameters  can  be  expressed 
through  the  use  of  small  parameters  £,  n,  C  as  follows: 

t^i e  )«j>»c  i  +iP(  0 )) 

e  >— l  +170(0 ))  (6) 

m. (  0  )— m,4C  1  +£A/(  0  )) 

where  P(e),  Q ( 0 )  and  M(6) — are  periodically  stable  random  processes 
with  zero  average  values.  We  defined  their  square  deviations  to  be 
1.  They  satisfy  the  following  periodic  conditions: 

Pi- *)-/»(*)»  Q(-  «)-Q(*) 

(y(- «)=.£?'(»), 


In  the  following,  we  used  the  periodicities  of  these  functions 
to  expand  these  functions  into  Fourier  series  using  the  spectral 
method:  «o 

p c«)»  2  /,«"* 


I  m  -99 


0(9)-  2 

/  *  -  oo 


(7) 


Mi  0  )—  2  h'e‘" 

I  X  -OO 


where  g.i—rft  h.t=h*,  are  the  random  Fourier  coefficients. 

*  represents  the  complex  conjugate. 
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Let  us  assume  that  the  random  processes  P ( 6 ) ,  Q ( 0 )  and  M(e)  are 
ergodic  processes  and  their  self-correlation  functions  can  be  obtained 
as  follows:  /?,(  t  )-£CP(  0 )P(  0  +  *  f  “  a  2  /«*"'  S 

I  ■  -  90  V  m  -90 

-  2  2  I/*!*®*  (8) 

I  9-99  I  9  l 

When  x  *  0,  the  self-correlation  function  is  equal  to  the  square  dev¬ 
iation.  Furthermore,  we  have  already  defined  the  square  deviation  of 
P(e)  to  be  1,  therefore, 

oo 

Rpi  0  )  —  2  2  1 

I  •  1 
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(8*) 


Similarly,  we  could  obtain  Rq(i)  and  Rm(-r). 

Let  us  assume  that  the  standard  deviations  of  the  random  pro¬ 
cesses  P(0),  Q(0)  and  M(0)  are  0  ,  oq  and  om>  respectively.  Then, 

from  the  above  analysis,  we  know  that  the  small  parameters  £,  n  and 
C  are  equal  to  these  standard  deviations,  respectively,  i.e., 

*  ,  (9) 

i  ’I  -0,1  t 

Prom  equations  (8)  and  (8'),  we  can  see  that  the  Fourier  coeffi¬ 
cients  such  as  | f , |  of  a  stable  ergodic  random  function  such  as  P(0) 
can  be  determined.  In  fact,  from  equation  (8*),  we  can  see  that: 

|/i|<o.707  (I  — J  ,2  ,-0  .  Furthermore,  with  regard  to  a  defined 

structure,  the  self-correlation  function  Rp(T)  can  be  experimentally 

obtained.  As  long  as  R  (t)  is  known,  then  the  numerical  value  of 

P 

| f ^ J  can  be  obtained  from  equation  (8). 

In  [4],  the  self-correlation  function  was  artificially  assumed 
to  have  the  following  form: 

/?,(*)-* "M  (10) 
when  the  parameters  a  could  be  determined  In  combination  with  the 
experimental  results.  From  the  known  equation  (10),  it  is  also  poss¬ 
ible  to  determine  the  actual  numerical  value  of  | f ^ ( . 

IV.  FREE  VIBRATION 


In  the  differential  equation  (5),  when  the  exciting  force  is 
zero  and  damping  is  neglected,  we  get 


-  WX'Y  + 


o 


The  boundary  conditions  are 


(11’) 


By  considering  the  randomness  of  the  structural  parameters,  q,  p 
and  mx,  and  by  substituting  equation  (6)  into  (11),  we  get 


-$JC(1  +  **?)*' 3' +C/>i(l  +tM))X-0 


The  obtained  equation  (12)  is  the  free  vibration  differential  equa¬ 
tion  of  a  periodic  random  parametric  system.  Because  the  variations 
of  the  random  parameters  are  infinitesimal  quantities  and  the  differ¬ 
ential  equation  contains  small  parameters  £,  n,  S,  it  is  possible  to 
find  its  solution  using  a  perturbation  method.  For  this  purpose,  the 
frequencies  and  vibration  modes  are  expressed  as  power  series  of  the 
small  parameters  £,  n,  and  5.  As  a  first  order  approximation,  let 
us  keep  the  first  order  terms  of  these  small  parameters: 

a>*  ■■  X*  +  4^*  +  ^v*  +  tP*  + . 

X(  9  )=  y  (  0  )  +  £«(  0  )  +  Tl«(  9  )  +  Cw<  9  )  + (13) 

By  substituting  them  into  equation  (12),  neglecting  higher  order  infin¬ 
itesimal  quantities,  and  making  the  coefficients  of  the  small  para¬ 
meters  £ ,  n ,  ?  zero,  we  get  the  following  equations: 


+  =  0 

-qlu'  u  -  -(PIP-*')  y 

-9iw''  +  (pi-A‘)  v  =  ql{Qy”+Q'y'  )  +  vly 
-  qlw’  +  (Pi - x1)  u>  -  ( X*A/ + p‘)  x 

The  first  equation  in  (1*0  is  a  well  defined  equation.  By  using 
the  boundary  conditions  to  solve  this  differential  equation,  we  obtain 
the  frequency  X2  as 

m  1  (15) 

where  m  is  an  integer  and  its  physical  meaning  is  the  nodal  diameter 

number  of  the  vibration  mode.  The  vibration  mode  y  (9)  is 

m 

+  (16) 

* 

where  A  and  A_  are  conjugating  coefficients  to  be  determined.  We 
m  m 

can  formulate  the  vibration  modes  in  such  a  way  so  that  these  coeffi¬ 
cients  can  be  determined.  Considering  the  following  orthogonal  con¬ 
ditions  : 


^  ym(  0  )  y»(  0  )^0“&« 


l  when  m  “  "  ) 
0  when  m  n  ) 


From  this  we  get 
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The  second  equation  (14)  contains  random  variables.  In  order 
to  solve  this  differential  equation,  we  expanded  the  function  u  (0) 
to  be  obtained  according  to  the  vibration  mode  y  (8)  by  letting 

u.<8)=  f|  6-,*(0>  (19) 

r  »  1 

Substituting  it  into  the  original  function  multiplying  y  (8) 
on  both  sides  of  the  equation,  and  integrating,  we  obtain  the  follow¬ 
ing  by  using  the  orthogonality  condition: 


+  Py-y*d 8 


when  m  =  n,  the  frequency  y2  obtained  from  equation  (20)  is 


y’  =  P«  J  =  (21) 


where  the  values  of  A  and  A*  can  be  given  from  equation  (18)  and 

m  m 

f2„  and  f*  from  equation  (7).  Let  us  assume 
mm 

Am-\AJe“~, 


then 


Pi  -  pil/i-|cos(  0,.  -  2a.) 


(21'  ) 


When  m  ^  n,  the  coefficient  b  „  of  the  vibration  mode  can  be 

mn 


obtained  from  equation  (20)  as 


Therefore, 


i  f l22' 

oo 

**-(  ® )—  2 


Following  a  similar  method  as  in  the  previous  section,  we  can 
obtain  solutions  for  the  third  and  fourth  equations  in  (14)  and  the 
frequency  obtained  is : 

v*  =  J  *  g  -ginr!p,.|cos(Y,.-2a.) 

(24) 

Pi--XiJ  *  n  Myic/0- -X2|A,Jcos(&,.-2a.) 

where  Am-\Am\'“-%  htm-\hlm\e‘'*-. 

The  expressions  of  the  vibration  modes  are: 


.  K* 
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(25) 


Vm(  9  )  “  ^  ®  ) 

n  *  l 
■  %  w 


w*(  9  )*=  S  £/--v-<9) 

»  «  1 

n\m  V 


Furthermore 


c»«”—  cMt  dmm*m—dm 


(26) 


V.  THE  NATURAL  FREQUENCY,  ITS  AVERAGE  AND  SQUARE  DEVIATION 


The  natural  frequencies  of  various  orders  of  the  system  can  be 
obtained  by  combining  equations  (13),  (15),  (21)  and  (2^).  The 
natural  frequency  of  mth  order  of  the  system  is: 

©i-xi+^l+nvi+CPi 

-(/tf+«‘d)  +ifil  J  *  ,  Pyi  dQ+  nrf  J  ‘  _  q/j </e-4XiJ  *  w  Myido  (27 ) 

-  (pj  +  «*9i)  +  |pil/,  Jcos(0,,  -  2a.)  -  T!m,gtl0t.|cos(  Y ,.  —  2P.) 

-  4  (pi  +  «*g#)|X,.|co8(&t.  — 2a.) 

In  searching  for  the  average  value  of  the  natural  frequency,  we 
have  to  consider  that  the  average  values  of  the  functions  P,  Q  and  M 
are  zero: 

£CF(9))-0»  £(Q(9))-0>  ECA/(  9  ))  =  0 


Hence,  the  average  value  of  the  natural  frequency  is: 

EOiD-E  CXi)  +  4/»jJ“.  ECF(9»  yW9+t?sj  EC(?(  9  ))yiV9 

-  tXi  J  *  #  ECAf(  9  ))  yldQ  -  Xi 


(28) 


The  average  square  value  of  the  natural  frequency  cu^  can  be  estimated 
in  the  following: 


E  C («0‘)  -  EC(Xi)*)  +  Vpi  J  *  „  J  "  „  E C  P  (9.)  P  (9,))yi(91);-i(9t)rf91rf9, 

+  T|‘  9,4  J  -  «  J  -  *  ^CO(®i)<?f9,)).v:i(0l).v:*(9lW0!,/9, 

+  t*Xi  J  *  ,  J  !  ,  E(Af(0,)Af(0t))yi(O,)yi(0t)d0|</0, 

+  24*1  Pi  <?i  J  "  ,  J  *  B  ECP(O,)O(9l)]jr£(Ol).v:t(Ol)</0l</Ol 

+  244  Pi  >.ij*  „  |  B  EC  F(0,)A/(0I))yi(0l)yi(O:)rf01d0t 
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+  21C  <jS  XiJ  *  _  ^  £,CO(01)A/(0l):>-  (O,)^i(O,)cf0,d0, 

+2>i  |  fcpij!,  £CP(  e  )Dyi(  e  )<fe+n<2ij* a  £C0(  e  9  Me 

r »  )  (29) 

-ttij  .  ,  £(M(  e  )Dyi(  9  M9  J 

The  square  deviation  of  the  natural  frequency  is: 

o’.  -E  C(®i)‘3 =  Y  {^il/,-l*  +  T»,m49ilfl«-l1  +  ttX:iftiJ‘ 

— 244plXil/» J  |At Jcoe(  P,«. — btm)  J 

+  2n4m‘qJXilfl,.!|Ai-|cos(Y„.- 5,.)} 

The  above  equation  indicates  that  th ?  standard  deviation  of  the 
natural  frequency  is  the  vector  sum  of  the  three  following  vectors 
(Figure  2): 

4pl/»-»  -tXlAt.  (31) 


The  standard  deviation  of  the  natural  frequency  can  be  estimated 
by  the  following  approximation: 

°.*<4/»JI/»J  +  WgilStJ  +  4X2|/i,-l  (32) 

We  can  see  that  the  standard  deviation  a  2  of  the  natural  fre- 

OJ 

quency  is  on  the  same  order  of  magnitude  as  the  small  infinitesimal 
quantities  £,  n»  4,  which  is  an  infinitesimal  quantity  on  the  same 
order  of  magnitude  as  the  standard  deviations  a^,  o  ,  am  of  local 
frequency  scatter. 


Fi|.  2  Standard  dtviatiou  of  satiral  frtqutncitt 
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VI.  NATURAL  VIBRATION  MODE  AND  ITS  PHASE  ANGLE 


The  vibration  modes  of  various  orders  of  the  system  can  be 
obtained  by  combining  equations  (13),  (16),  (23)  and  (25).  The  mth 
order  natural  vibration  mode  is 


Xm(  e  )-y.(  0  )  +  l«-(  0  )  +  *,„.(  0  )  +  tw-(  0  ) 

•  (33) 

-y-(0)+  2  C-+W-W0) 

*  -  1 
ii 


Due  to  the  complexity  of  the  vibration  mode  expression  (33), 
we  can  see  that  the  nodal  diameter  of  the  vibration  mode  is  changing 
from  a  symmetric  distribution  to  an  asymmetric  distribution  (Figure 

3). 


It  can  be  proven  that  the  vibration  mode  (33)  satisfies  the 
orthogonality  condition: 


j".,  XmXm'  ' 


{lwhen  m  =»m' ) 
Owhen  m  ) 


(34) 


As  for  the  phase  angle am  of  the  natural  vibration  mode,  when  the 
structural  parameters  are  constant,  the  phase  angle  am  is  a  value  to 
be  determined  by  the  Initial  condition.  When  the  structural  para¬ 
meters  are  not  constants  but  random  quantities,  the  natural  vibration 
mode  phase  angle  am  is  determined  by  the  variations  of  the  structural 
parameters  P(0),  Q(6)  and  M(9)  and  Is  not  related  to  the  initial  con¬ 
dition. 

In  fact,  in  solving  for  the  latter  three  equations  in  (14),  the 
right  hand  side  of  the  equation  should  not  contain  a  "long  term" 
term  with  spatial  frequency  m.  Otherwise,  a  solution  which  Increases 
with  time  will  be  found  which  is  not  rational. 

In  order  to  make  sure  that  the  solution  Xm(  0  )-y-+4«.+Tit>*+Cw. 
does  not  contain  any  term  which  increases  with  time,  we  must  let  the 
coefficients  of  the  following  equation  be  zero: 
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Fi*.  3  Nodal  di»®*Uf*  of  natural  mod** 
<«-  3) 

- - »f»m«  tricali  - **ym®  e  t  ricul . 


Fij.  4  Campbell  diHt*“ 
q - violent  •— e.k  retonanca. 


{--Mr  l/t  Jsin  (0*.  —  2a.)  +  Tlm,|ff,Jsin( Yt„  —  2«„) 

+CAI  |  A,J  sin  ( btm  —  2a»  )  j  sin  (  m  0  +0“  0 

After  expansion,  we  can  obtain  the  formula  to  calculate  the 
phase  angle  of  the  vibration  mode: 

„  - 1  flil/.Jsin0..+ Wqjlg,  Jstnr,.  ±  t  (35) 

tg2a.  -  —  ^JI/l.lcosP,.  +  ,lm*qJ|g,.|cosY +  tAil  A.  Jco9&,» 

VII.  FORCED  VIBRATION 


In  equation  (1),  let  us  assume  that  vtie  exciting  force  f(9,t) 
contains  many  harmonics: 


/(9,o-  S  /. c e .  o-  S  W-F*( 9 )e" 


Now,  let  us  only  consider  the  k  harmonic  of  the  exciting  force, 
i.e.,  on  the  right  hand  side  of  equation  (1)  we  only  take 

/,(«.  I  )-m.,F,(9)ei‘Q'  (37) 

For  the  convenience  of  writing  without  losing  the  generality, 
let  us  assume  that  the  spatial  distribution  of  the  exciting  force 

f*  h 

and  the  k  harmonic  of  the  vibration  mode  have  the  same  phase  angle, 
which  is  to  assume  that 

F*(9)-F,y,(0)  (38) 


,  a 

,  .  .  - 

'v- 
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Let  us  assume  the  special  solution  to  the  differential  equation 
(1),  i.e.,  the  forced  vibration  is 

*(0,  i  )-.Y(  9  V*°'  (39) 

Substituting  it  into  equation  (1),  we  obtain  the  differential 
equation  of  forced  vibration  as 

-(9** ')'+[/»*-  (Ml)*  ].Y+«efcQ*-F»v*(0)  (40) 

In  the  following,  in  order  to  simplify  the  derivation,  let  us 
assume  mx(0)  *  mxQ  *  const  which  is  to  neglect  the  nonuniformity  of 
mass  distribution  by  assuming  £M(0)  =  o. 

Let  us  assume  that  the  solution  to  the  differential  equation 
(40)  can  be  expressed  as  the  sum  of  the  various  orders  of  natural 
vibration  modes: 

*<•>-  2  *-*-<•> 


where  the  natural  mode  X  (9)  can  be  found  in  equation  (33).  Substi- 

m 

tute  equation  (41)  into  the  differential  equation  (40).  Multiply 
both  sides  by  Xm(0)  and  then  integrate.  By  using  the  orthogonality 
condition  (17)  and  neglecting  higher  orders  of  the  infinitesimal 
amounts,  we  get 

(  F*  (when  "» *  fc  ) 

<«i-(  *Qy+UkO)a.-< 

\  F,(^..,  +  ^.*)(whenm^fc)  (42) 


(  F,  (vd 

I  F»(  £6„,»  +  tc«»)  1 


VIII.  RESONANCE  AMPLITUDE  AND  ITS  AVERAGE  AND  SQUARE  DEVIATION 

The  following  is  a  consideration  of  the  resonance  amplitude. 
Please  refer  to  the  Campbell  diagram  in  Figure  4. 

(1)  when  the  "triple  point"  condition  is  satisfied,  i.e.,  when 
and  m  *  k,  from  equation  (42): 

itkQ 

This  is  the  quantity  of  a  defined  mode.  Under  small  damping  condi¬ 
tions,  the  amplitude  is  very  large  and  a  very  intense  resonance  is 
obtained . 


(2)  when  the  "triple  point"  condition  is  not  satisfied,  i.e., 
when  ,  but  m  ^  k.  From  equation  (k2): 

w-ew(^+Tk-}  (l424) 

At  this  time,  the  average  of  the  vibration  amplitude  is  zero: 


S'W’-sfe-  ))/.>•: }./o-  o  („5) 


The  square  deviation  of  the  resonance  amplitude  is 

J-*  fc/,(o,)Q(eI)iy-(eI)y.ce,)y-(o,)y«(e,)i/e1de, 

+  1^1' b  |*n  £CO(e,)O(e.)iy:(e1)^*(9,)y:(0t)^(fli)Jfl.‘/0« 

+an- J { "  „  J " , £,CP(0.)Q(0«)3y-(0.)y*(0.)y-(0*)y*^*)£,0>d0‘}  (46 } 

The  calculated  results  are: 

a‘ "  2  (n,-V),(«W)1  !*‘  +  T,*m< 


~2tW  +  V-4-|/..,|« 

9i  <?• 

+  i>W|*..,|*  +  24W  -g-  \fm.t\\gm.>\  oo8(p..,-Y..,)  } 
0*  I 


(47) 


From  the  above  complicated  expression,  we  can  see  that  the 

ratio  of  the  standard  deviation  a  of  the  resonance  amplitude  and  the 

a 

resonance  amplitude  |amoj  under  the  triple  point  condition  can  be 
considered  as  the  vector  sum  of  the  following  vectors  (Figure  5): 

4  /•»«!  4  ~^hf m.»t  (48) 

This  can  be  estimated  approximately  using  the  following  equation: 

4 + +^4 -^-1/— »I+Wls..»i^  |  (^9) 


From  this  we  can  3ee  that  this  ratio  is  on  the  same  order  of 
magnitude  as  the  small  parameters  5,  q  and  the  standard  deviation  op, 
aq  of  the  scatter  of  the  local  frequencies. 
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Fif.  S  Squrt  dtrittioas  of 


refoniace  amplitude 


IX.  EXAMPLE 

With  respect  to  the  mechanical  model  shown  in  Figure  1,  let  us 
choose  the  local  frequency  as: 

A-110X2«,  ft  *■  22  x  2  * 

Usually,  a  turbine  machinery  factory  establishes  that  the  fre¬ 
quency  scatter  of  the  blades  on  a  turbine  cannot  exceed  ±U%.  Now, 
let  us  take  1/3  of  this  value  as  the  standard  deviation  of  the  fre¬ 
quency  scatter,  i.e., 

t  -  TJ  -a,-ot--^yi-=0. 01333 

Assuming  that  the  self-correlation  function  of  the  random  struc¬ 
tural  parameters  has  the  form  of  that  of  equation  (10),  let  us  also 
choose  a  *  1.  Using  equation  (49),  let  us  calculate  the  extremes  of 
the  ratio  of  the  standard  deviation  of  resonance  amplitude  and 
the  resonance  amplitude  jamQ|  under  the  triple  point  conditions  at 
various  orders  of  harmonics  k  of  the  exciting  force,  and  different 
main  vibration  mode  nodal  diameter  numbers  m  are  shown  in  Table  1. 

From  Table  1,  we  can  see  that,  under  the  condition  of  the  main 
mode  with  no  nodal  diameter  (m  *  0)  and  the  first  order  harmonic 
exciting  force  (k  *  1),  oj  |o«,|<13.58*  Under  other  conditions, 


Table  1  Limit  values  of  ratio  o./|o«,| 


harmonic  nodal  dia 
number  meter 

m  •  0 

no. 

M  >  1 

*  «  2  | 

*"* 

•  •  4 

k  *  0 

1 

0.1413 

00239 

0.00*2 

0.0046 

*  •  1 

0  135* 

i 

0-0467 

0.0112 

0.0055 

k  -  2 

0  020* 

00365 

I 

0.0271 

0.00*5 

*  -  1 

0.0067 

0  0086 

0.0231 

1 

0.022* 

k  •  4 

0  002* 

0  0035 

0  0060 

0.01*0 

1 

*  -  5 

0.0015 

0  001* 

0  0026 

0.0051 

0.0177 

0  000* 

0.0010 

0.0013 

0.0023 

0.004* 

this  ratio  is  extremely  small  which  can  be  neglected.  References 
[6-7]  have  reported  the  experimental  results  on  the  machine.  They 
believed  that  the  resonance  amplitude  of  a  circumferentially  con¬ 
nected  blade  group  was  approximately  ~j  of  the  resonance  ampli¬ 
tude  of  free  blade.  This  result  is  on  the  same  order  of  magnitude 
as  the  calculated  result  based  on  the  above  model. 

X.  CONCLUSIONS 

1.  For  a  periodic  dynamic  system  with  random  structural  para¬ 
meters,  to  use  a  spectral  method  to  find  the  solution  is  a  feasible 
and  convenient  method.  The  random  functions  of  the  structural  para¬ 
meters  are  expanded  into  Fourier  series  to  facilitate  the  calculation 
of  the  natural  frequency,  natural  mode  and  resonance  amplitude  of  the 
system,  and  to  estimate  their  average  values  and  the  square  deviations. 

2.  The  analytical  calculation  with  regard  to  a  periodic  dynamic 
system  with  random  structural  parameters  indicated  that  its  vibration 
characteristics  have  some  special  features  different  from  those  of  a 
uniform  periodic  structure. 

(1)  the  natural  frequencies  of  various  orders  of  a  uniform  struc¬ 
ture  are  fixed  quantities,  while  those  of  a  random  structure  are  ran¬ 
dom  quantities.  The  standard  deviation  of  the  mth  order  natural  fre¬ 
quency  is  only  related  to  the  2mtil  order  Fourier  coefficient  of  the 
random  parameter.  Furthermore,  it  is  the  vector  sum  of  the  standard 
deviations  of  several  deviations. 


(2)  the  natural  vibration  modes  of  a  uniform  structure  are  of 
harmonic  waveforms  and  all  the  nodal  diameters  are  uniformly  distri¬ 
buted.  For  a  random  structure,  in  addition  to  the  harmonic  wave¬ 
forms  of  the  natural  frequencies  of  the  various  orders,  there  are 
various  orders  of  harmonic  waves  with  complicated  shapes  and  uneven 
nodal  diameter  distribution.  However,  all  the  natural  vibration  modes 
satisfy  the  orthogonality  condition. 

(3)  the  phase  angles  of  the  natural  modes  of  a  uniform  structure 
are  parameters  to  be  determined  by  the  initial  condition.  However, 
the  phase  angles  of  the  vibration  modes  of  a  random  structure  are 
related  to  the  respective  structural  parameters  which  are  not  related 
to  the  initial  condition. 

(4)  for  a  uniform  structure,  the  "triple  point"  condition  must 
be  satisfied  to  create  resonance.  Under  the  condition  that  the 
"triple  point"  condition  is  not  satisfied,  it  is  impossible  to  create 
resonance.  For  a  random  structure,  when  the  "triple  point"  condition 
is  satisfied,  very  strong  resonance  will  be  created.  However,  even 
when  the  "triple  point"  condition  is  not  satisfied,  as  long  as  &m-kQ, 
even  though  k  /  m,  it  will  create  a  weak  resonance.  In  this  paper  we 
estimated  the  square  deviation  of  the  weak  resonance  amplitude.  This 
square  deviation  is  related  to  certain  orders  of  the  Fourier  coeffi¬ 
cients  of  the  structural  parameters. 


96 


REFERENCES 


[1]  Huang  Wenhu,  Zhao  Yuchang,  Den  Liancbao.  "The  Analysis  of 
Vibration  of  a  Circumferentially  connected  Blade  Group". 

Harbin  Turbine  Factory,  Technical  Promotion  Office,  1974. 

[2]  Huang  Wen  bu,  Deng  Lianchao,  Zhao  Yuchang.  "The  analysis  of 
Vibration  of  one  Kind  of  Periodic  Structures".  Journal  of 
Harbin  Institute  of  Technology,  vol.  3,  1979,  PP  11-30. 

[33  Huang  Wenhu.  Free  and  Forced  Vibrations  of  Closely  Coupled 

Turbomachinery  Blades,  AIAA  Journal,  vol.  19,  no.  7,  July  1981 

[4]  Hoshiya,  M.  and  Shah,  H.  C.  Free  Vibration  of  Stochastic  Beam' 
Column,  EM.  Div. ,  ASCE,  Aug.  1971. 

[5]  Sing  Kushen.  "Random  Vibration  Analysis".  Earth  Publications 
1977. 

[6]  Yampol ’skaya,  R.G.,  and  Arkad’yev,  D . A .  Determination  of  the 
vibration  Characteristics  of  the  Blading  with  Damper 
Connections.  Energomashchinostroyeniye ,  No.  11,  1965. 

[7]  Arkad'yev,  D.A.  Influence  of  Couplings  on  the  Vibration 
Strength  of  Turbine  Blades.  Energomashchinostroyeniye, 

No.  3,  1968. 


97 


A  SPECTRAL  APPROACH  FOR  ANALYZING  THE 
VIBRATION  OF  A  PERIODIC  STRUCTURE  WITH 
RANDOM  PARAMETERS 


Huang  Wenhu 

{Harbin  Institute  of  Technology ) 

Abstract 

In  •  periodic  structural  system  such  as  blades  in  a  circumferentially  clo¬ 
sed  packet  on  a  disk  of  turbo-machinery*  the  natural  frequencies  of  individu¬ 
al  blades  can  be  randomly  different  from  one  another.  Prom  this  arises  the  pro¬ 
blem  of  vibration  analysis  of  a  periodic  structure  with  random  parameters. 
There  is  lack  of  general  method  for  solving  the  differential  equations  with  ran¬ 
dom  parameters.  This  paper  describe  ..  ;  ectral  approach  for  analyzing  the  vi¬ 

bration  of  a  periodic  structure  with  random  parameters.  Suppose  the  standard 
deviations  of  random  structural  parameters  are  small  so  that  a  perturbation  me¬ 
thod  can  be  used  to  reduce  the  differential  equation  with  several  random  pa¬ 
rameters  to  several  differential  equations  with  one  parameter,  and  then  these 
differential  equations  may  be  solved  one  by  one.  Suppose  the  spatial  distribu¬ 
tions  of  the  random  structural  parameters  are  ergodic.  and  for  concrete  stru¬ 
cture  these  distribution  functions  and  their  correlation  functions  can  be  deter¬ 
mined  by  experiments.  It  is  suggested  in  this  paper  to  expand  these  spatial  dis¬ 
tribution  functions  of  random  parameters  into  Fourier  Series.  Ami  then  the  re¬ 
lation  between  these  Fourier  coefficients  and  the  correlation  functions  is  esta¬ 
blished  so  that  these  Fourier  coefficients  can  be  detei mined  by  several  ways. 

In  this  situation,  these  differential  equations  with  random  parameters  can  be 
solved.  Thus  natural  frequencies  of  the  structure  are  then  obtained,  and  their 
standard  deviations  are  estimated.  Also,  the  expressions  of  natural  modes  are 
given,  the  orthogonality  of  natural  modes  is  proved,  and  it  is  shown  tl  at  the 
phase  angles  of  natural  modes  are  not  arbitrary.  Finally  the  special  conditions 
of  resonance  of  periodic  structure  with  random  parameters  are  discussed.  It  is 
shown  that  a  violent  resonance  occurs  when  the  number  of  harmonic  of  excit¬ 
ing  force  is  equal  to  the  number  of  nodal  diameters  of  natural  modes,  and  only 
a  weak  resonance  appears  when  these  two  numbers  are  not  equal.  This  pheno¬ 
menon  does  not  exist  in  the  case  of  structures  with  homogeneous  parameters. 
The  standard  deviations  of  amplitudes  of  weak  resonance  are  estimated.  Nu¬ 
merical  examples  show  that  the  calculated  results  have  the  same  order  as  the 
experimental  results  in  literature. 


A  NEW  METHOD  OP  FINITE  ELEMENT  STRUCTURE  DISCRETIZATION 


Liang  Guowei* 

Northwescern  Polytechnical  University 


ABSTRACT 

In  this  paper,  the  iso-parametric  element  geometric 
interpolation  method  was  used  to  automatically  generate 
the  nodal  coordinates  of  all  the  elements  based  on  a  small 
amount  of  input  data.  A  "chessboard"  mesh  was  used  to 
facilitate  the  numbering  of  elements  and  nodal  points. 
Simultaneously,  a  "front  solver"  method  was  used  to  solve 
the  equations  to  simplify  the  numbering  program. 

This  method  has  the  advantages  of  little  input  data, 
ease  of  changing  the  mesh  and  the  sufficient  accuracy  of 
the  boundary  nodal  coordinates.  It  has  been  used  for  the 
two  dimensional  mesh  of  an  axial  symmetric  body  and  the 
three  dimensional  mesh  of  a  turbine  blade.  The  results 
were  satisfactory. 


I .  INTRODUCTION 


In  the  finite  element  calculation  work,  the  work  load  for  the 
preparation  of  input  data  is  very  large.  It  mainly  Involves  the  deter¬ 
mination  of  nodal  point  coordinates  and  number.  If  these  data  are 
obtained  manually,  it  Is  very  easy  to  make  mistakes.  Therefore,  this 
will  cause  difficulties  in  testing  the  program  and  wasting  computer 
time.  Hence,  automatic  discretization  becomes  an  actual  problem  to 
be  resolved  in  the  finite  element  computation  work. 

In  the  references  [1-4]  both  here  and  abroad  in  the  70's,  various 
methods  to  automatically  form  meshes  were  discussed;  which  were 
limited  to  planar  triangular  elements. 

As  we  all  know,  under  the  same  nodal  point  number  the  accuracy 
of  a  rectangular  element  (or  hexagonal  element)  Is  higher  than  that 
of  a  triangular  or  (tetrahedral)  element.  Therefore,  it  Is  very 
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imperative  to  develop  a  rectangular  iso-parametric  element  which 
is  applicable  to  the  automatic  generation  of  two-dimensional  and 
three  dimensional  meshes.  Furthermore,  in  the  aforementioned  litera¬ 
ture,  the  emphasis  was  focused  on  the  automation  of  the  division 
of  the  element  and  the  formation  of  the  nodal  point  coordinates.  As 
for  the  automatic  numbering  of  the  nodal  points  and  elements,  it  was 
seldom  considered  in  connection  with  the  storage  problem.  In  some 
references  [5-6],  although  some  methods  for  renumbering  the  nodal 
points  and  reducing  the  bandwidth  were  proposed,  yet  such  methods 
would  complicate  the  program. 

In  1978,  the  author  proposed  a  method  to  automatically  form  the 
nodal  point  coordinate  using  an  interpolation  function  [7]  together 
with  a  "front  solver"  method  to  solve  the  linear  equations  by  consi¬ 
dering  the  convenience  of  automatic  numbering  of  the  nodal  points. 

The  numbers  of  the  nodal  points  were  allowed  to  be  discontinuous 
to  simplify  the  program  as  well  as  to  save  the  internal  storage. 

During  the  numbering  of  the  elements,  some  programming  techniques  were 
used  to  minimize  the  "wave  front".  After  the  discretization  of  the 
mesh,  a  huge  amount  of  data  was  required  to  be  checked.  Undoubtedly, 
it  is  very  time  consuming.  This  method  displays  these  data  using  the 
method  most  convenient  for  checking  which  consequently  saves  a  lot 
of  time. 

II.  MAIN  POINTS  OF  THIS  METHOD. 

1 .  Automatic  formation  of  nodal  point  coordinates 
using  the  interpolation  function 

Usually,  geometric  interpolation  is  used  in  an  iso-parametric 
element  to  transform  a  rectangular  (or  hexagonal)  element  into  a  curve 
(curved  surface)  element.  The  coordinates  of  an  arbitrary  point  P 
of  the  curve  (or  curved  surface)  element  can  be  obtained  by  the  nodal 
coordinate  Interpolations  of  the  elements: 


*  =  2  V/U,  n,  t ;  •  *, 
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|  (A) 

y  ~  J  -V/(  4 ,  m  .  t )  •  >v 

» 

*  =  £-V,(  4.  ’1.  t )  •  zt 


where  (x.^,  yi5  zi)  is  the  coordinate  of  nodal  point  i,  (£,n,s)  is 
the  local  coordinate  of  the  point  P,  and  N,(4,n,4)  is  the  shape 
function  corresponding  to  the  nodal  point  i  as  shown  in  Figure  1. 

Through  the  transformation,  the  original  linear  mesh  (the  iso  £, 
iso  n  lines  on  the  left  of  Figure  1)  is  transformed  into  the  curve 
mesh  (the  iso  £  and  iso  n  lines  in  the  x-y  coordinate  system  on  the 
right  of  Figure  1).  If  this  iso-parametric  element  is  the  component 
we  want  to  analyze  (or  part  of  a  component),  then  the  large  amount 
of  nodal  point  coordinates  in  the  mesh  can  be  obtained  by  interpolat¬ 
ing  the  small  amount  of  original  nodal  coordinates  (x,,  y,  ,  z^)  using 
equation  (A). 


Fi*-  1  2-tad  3-dimeatioatl  ao-ptrtat)rit  eleaeat* 
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2 .  Actual  procedures  of  nodal  point  coordinate  Interpolation 


(1)  Based  on  the  boundary  shape  of  the  component,  it  is  divided 
into  several  zones.  For  each  zone,  a  proper  parent  element  for  each 
zone  was  chosen; 

(2)  in  the  body  coordinate  system,  based  on  the  boundary  shape 
of  the  component,  the  nodal  point  coordinates  (x.,,  y^,  z.)  of  the 
parent  element  was  chosen; 

(3)  based  on  the  degree  of  closeness  of  the  discretized  mesh,  the 
local  coordinates (£, n ,? )  of  any  nodal  point  on  the  parent  element  in 
the  mesh  in  that  zone  are  determined; 

(4)  from  equation  A),  the  coordinates  (x,y,z)  of  all  the  nodal 
points  of  the  mesh  in  that  zone  can  be  interpolated. 

Now,  let  us  use  the  disk  in  Figure  2  as  an  example  to  explain 
the  above  procedures. 

First,  based  on  the  shape  of  the  disk,  let  us  divide  it  into  18 
zones.  The  principle  of  zoning  is  that:  a  different  boundary  curve 
should  be  divided  into  two  zones.  For  example,  the  a-b-c  boundary 
in  the  figure,  a-b  is  a  line,  b-c  is  a  circular  arc,  then  a-b-c  6 

should  be  divided  at  point  b,  and  a-b  is  in  zone  7  and  b-c  is  in  zone 

8. 


For  each  zone,  a  parent  element  is  selected  and  its  nodal  points 
are  determined  by  the  boundary  shape.  For  zone  7,  the  boundaries  are 
straight  lines  and  a  four  nodal  point  parent  element  is  chosen  (no.  1 
in  the  figure)  in  order  to  ensure  the  accuracy  of  boundary  interpola¬ 
tion.  For  zone  8,  it  is  necessary  to  choose  the  no.  2  parent  element 
to  ensure  the  sufficient  accuracy  of  the  interpolated  b-c  boundary. 

The  principle  of  selecting  the  parent  element  is  that:  under  the 
premise  of  assurance  of  sufficient  accuracy  of  boundary  interpolation, 
the  nodal  points  should  be  as  few  as  possible  in  order  to  reduce  the 
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input  coordinate  data.  In  the  mean¬ 
time,  it  should  be  ensured  that  the 
mesh  in  the  neighboring  zones  is  conti¬ 
nuous  . 

Next,  on  the  cross-section  of  the 
disk,  the  nodal  point  coordinates  (R^, 
z^)  in  this  zone  are  determined  based 
on  the  nodal  point  distribution  of  the 
selected  parent  elements  in  each  zone. 

For  corner  nodal  points,  their  coordi¬ 
nates  are  fixed.  As  for  nodal  points 
located  on  the  side,  there  is  room  for 
selection.  The  numbers  and  positions 
of  the  middle  line  nodal  points  would 
affect  the  error  of  the  interpolated 
boundary  curve  greatly.  Based  on  the 
experience  of  the  author,  for  a  second 
order  curve,  the  use  of  1-2  middle 
nodal  points  Is  sufficient.  Further¬ 
more,  a  uniform  location  distribution 
might  be  optimal.  For  a  complicated 
shape  such  as  the  turbine  blade,  two 
middle  nodal  points  are  used.  The  maximum  error  of  the  interpolated 
boundary  nodal  point  coordinates  is  within  ±0.15  mm.  With  regard  to 
stress  analysis,  it  can  be  considered  satisfactory.  In  this  case,  the 
middle  nodal  point  positions  may  not  be  uniformly  distributed.  Inci¬ 
dentally,  by  changing  the  position  of  the  middle  nodal  point.  It  is 
possible  to  change  the  density  distribution  of  the  mesh  which  is  one 
of  the  advantages  of  the  method.  However,  it  should  be  noted  that. 

If  the  distance  between  the  middle  nodal  point  and  the  corner  nodal 
point  is  less  than  1/4  (1  is  the  length  of  the  side),  then  the  mesh 
is  "wrinkled”  [8-9]. 

The  local  coordinates  (£,n,C)  of  any  nodal  point  in  the  mesh  are 
easily  obtained  in  the  parent  element  because  the  mesh  Is  uniform  in 
the  parent  element.  As  long  as  the  longitudinal  and  transverse  lines 
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zones  and  their  parent  dements 
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of  the  mesh  are  chosen.  It  is  not  difficult  to  determine  the  partial 
coordinates  of  any  nodal  point. 

Due  to  the  fact  that  the  distribution  of  the  nodal  points  of 
the  parent  elements  in  each  zone  is  irregular,  the  shape  function 
can  be  directly  derived  using  the  method  in  [10]. 

3 .  The  numbering  of  elements  and  nodal  points 

In  designing  the  program,  it  is  more  complicated  to  number  the 
nodal  points.  If  optimized  numbering  is  used  to  reduce  the  band¬ 
width,  it  would  complicate  the  program  further.  The  author  adopted 
the  "front  solver"  method  to  solve  the  set  of  linear  operations.  It 
has  two  advantages:  "front  solver"  method  itself  can  save  the  inter¬ 
nal  storage,  besides,  "front  solver"  method  does  not  have  a  bandwidth 
requirement  for  the  numbering  of  the  nodal  points.  Even  when  the 
nodal  point  sequence  is  not  continuous,  it  does  not  matter  at  all  [11]. 
This  makes  It  difficult  to  compile  the  program. 

The  procedures  to  number  the  nodal  points  are  as  follows: 
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COMPUTER  RESEARCH  ON  DYNAMIC 
CHARACTERISTICS  OF 
A  SYNCHRO  GENERATOR 

Qian  Zhenxicmg  and  Xu  Qiaobao 

( Beijing  Institute  of  Aeronautics  and  AstrOnau.'ics) 

Abstract 

A  P-transform  matrix  is  applied  to  transforming  the  time-varying  dyna¬ 
mic  equations  of  a  synchro  generator  into  so-called  P  Park's  equations, in  re¬ 
sult  the  state  equations  of  one  unit  generator  system  are  established.  Owing  to 
the  nonlinear  effect  of  the  voltage  regulator  and  the  lagging  of  the  a.  c.  exciter 
in  an  aircraft  brushless  a.  c.  generator,  the  forcing  function  of  the  main  synchro 
generator  C’,(  t  )  can  be  reduced  to  an  exponential  curve  with  a  lagging  ts  and 
an  equivalent  time  constant  T 

On  the  basis  of  the  state  equations  above  mentioned  eight  dynamic  chara¬ 
cteristics  of  the  synchro  generator  are  printed  out  separately  on  an  electronic 
digital  computer  by  four  different  methods,  analysis!  fourth-order  Rungo-Ku- 
tta  algorithm!  exponential  matrix  (*AT)%  and  network  topology.  The  computa¬ 
tion  accuracy  and  the  stability  region  of  th"S“  four  methods  are  analyzed  and 
their  applicable  range  is  established  in  conformity  with  their  advantages.  Then 

the  following  conclusions  are  drawn, 

1.  The  characteristic  roots  of  the  coefficient  matrix  at  no-load  sudden  short 
circuit  are  the  dynamic  parameters  of  the  synchro  generator,  i.e.  T»,  TJ.  'I  Tj. 

2.  In  symmetrical  operation  the  synchro  generator  can  be  simplified  as  an  one 
-order  inertial  link.. 

3.  The  dynamic  characteristics  of  the  a.  c.  exciter  (u,  T„)  have  rather  little 
effect  on  the  maximum  surge  current  of  the  synchro  generator  at  sudden  short 
circuit. 

4.  In  consideration  of  the  non-linearity  of  the  magnetic  circuit  in  the  synchro 
generator,  it  is  verified  that  the  inductor  flux  it  and  the  capacitor  charges  qc 
are  more  suitable  than  ii.  and  uc  to  be  taken  as  the  state  variables. 
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TECHNICAL  EXCHANGE  MEETING  ON  RELAYS  AND  CONTACTORS  OF  THE 
AERONAUTICAL  SOCIETY  OF  CHINA 

The  China  Society  of  Aeronautics  and  Astronautics  held  a  meeting 
entitled  ‘'Technical  Exchange  Meeting  on  Relays  and  Contactors”  in 
Zunyi  between  March  12-17,  1982.  The  meeting  was  organized  by  the 
315th  Factory  of  the  Third  Machinery  Department.  The  delegates  attend¬ 
ing  the  meeting  included  73  people  from  35  organizations  from  the  Third 
Machinery  Department,  the  Fourth  Machinery  Department,  the  Seventh 
Machinery  Department,  the  First  Machinery  Department,  the  Air  Force 
and  the  Navy. 

The  conference  received  over  30  papers  and  technical  reports. 

17  were  presented  in  the  conference  and  over  10  were  read  in  group 
meetings.  The  papers  and  reports  covered  a  wide  range  of  areas. 
Primarily,  they  are  the  results  in  theoretical  design,  new  product 
development  and  exploratory  directions  in  the  relay  field,  such  as  the 
development  of  the  computer  assisted  design  and  verification  of  the 
dynamic  parameters  of  a  magnetic  relay  system,  the  complete  verifica¬ 
tion  of  the  JKM  sealed  relay  series,  the  10  A  4  circuit  magnetic  main¬ 
tenance  relay,  the  0.5  A  2  circuit  miniature  magnetically  maintained 
relay,  the  square  TO-5  transistor  tube  case  sealed  relay,  and  the 
spherical  TO-5  relay  as  well  as  the  development  of  new  silver  oxide 
contact  materials. 

During  the  meeting,  the  delegates  conducted  a  panel  discussion  on 
the  problem  of  the  developmental  direction  of  the  relay. .  It  was  unan¬ 
imously  agreed  upon  that  aeronautical  relays  and  contactors  should 
develop  in  the  following  three  areas:  (1)  the  maneuverability  of 
attack  aircraft,  (2)  the  modernization  of  aircraft  control  system  and 
engine  control,  and  (3),  the  development  of  new  aircraft  power  sources. 
In  summary,  efforts  will  be  continued  based  on  high  reliability,  high 
sensitivity,  high  velocity,  low  load  capability,  low  power  assumption 
and  miniaturization. 

Finally,  the  meeting  suggested  that  the  next  technical  exchange 
meeting  be  held  in  1984. 
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ABSTRACT 

The  computer  aided  design  technique  is  an  important 
development  in  computer  applications  and  it  is  an  import¬ 
ant  component  of  computer  science.  The  special  language 
for  electronic  circuit  analysis  is  the  foundation  of  com¬ 
puter  aided  design  or  computer  aided  circuit  analysis 
(abbreviated  as  CACD  and  CACA)  of  simulated  circuits. 

Electronic  circuit  analysis  language  (ECAL)  is  a  com¬ 
paratively  simple  and  easy  to  use  circuit  analysis  special 
language  which  uses  the  Fortran  language  to  carry  out  the 
explanatory  executions.  It  is  capable  of  conducting  dc 
analysis,  ac  analysis,  and  transient  analysis  of  a  circuit. 
Furthermore,  the  results  of  the  dc  analysis  can  be  used 
directly  as  the  initial  conditions  for  the  ac  and  transient 
analyses . 

The  ECAL  language  describes  the  circuit  by  using  input 
statements  which  are  familiar  to  electronic  circuit  engin¬ 
eers.  As  for  the  allowable  elements,  in  addition  to  the 
regular  linear  elements,  such  as  resistors,  capacitors, 
inductors,  current  sources  and  voltage  sources,  nonlinear 
elements,  such  as  diodes,  transistor  triodes  and  nonlinear 
resistors,  are  also  permitted.  Hence,  it  is  capable  of  cir¬ 
cuit  analysis  for  both  linear  and  nonlinear  circuits. 

The  ECAL  language  uses  very  simple  output  statements  to 
control  the  output  form  of  the  resultant  analyzed  data.  It 
may  be  tables  or  figures,  depending  on  the  various  needs. 
Therefore,  it  is  a  very  useful  analytical  tool  for  engineers. 


INTRODUCTION 


Computer  aided  design  is  an  Important  aspect  of  computer  applica¬ 
tions  which  is  also  an  important  milestone  in  the  development  process 
of  computer  applications  from  elementary  to  advanced  stages.  As  a 
branch  of  computer  science,  the  computer  aided  design  technique  is  a 
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special  technique  which  is  the  simple  application  of  the  computer  to 
a  process  or  part  of  the  design  work.  Even  for  a  product,  in  spite 
of  the  fact  that  its  design  has  been  completed  by  a  computer  from  the 
beginning  to  the  end,  if  the  program  used  was  compiled  specifically 
for  the  design  of  this  product,  it  is  still  not  possible  to  claim 
that  the  computer-aided  design  technique  was  adopted.  In  other  words, 
computer  aided-design  technique  has  certain  conditions  as  its  label. 
These  conditions  are:  1)  there  is  a  suitable  special  language  for 
this  type  of  problem;  2)  there  is  a  data  bank  with  various  computa¬ 
tional  methods  and  applications  data;  3)  there  is  a  combined  software- 
hardware  system  with  a  certain  dialog  capability  between  man  and 
machine  to  allow  the  designer  to  interfere  with  the  work  of  the  com¬ 
puter  at  any  time.  This  allows  the  man  and  the  machine  to  do  their 
best  to  complete  this  task  fast  and  well. 

The  electronic  circuit  analysis  language  ECAL  was  developed  as  an 
important  component  of  the  computer-aided  design  technique  of  electron¬ 
ic  circuits.  By  taking  into  account  the  convenience  of  use  for  the 
circuit  designer  and  the  feasibility  of  a  man-machine  dialog,  its 
input  statement  form  is  simple,  which  is  consistent  with  the  custom 
of  the  circuit  analyst.  Therefore,  it  is  easy  to  learn  and  convenient 
to  use. 

The  electronic  circuit  analysis  language  ECAL  is  a  comparatively 
simple  and  easy  to  use  language  for  circuit  analysis  which  uses  the 
Fortran  language  to  carry  out  the  explanatory  executions.  It  includes 
a  set  of  statements  with  simple  structures.  The  statements  have  the 
capability  to  describe  the  circuits.  They  also  have  the  capabilities 
to  designate  the  analysis  range  and  the  output  form.  In  addition, 
they  are  also  able  to  control  whether  the  analysis  should  be  repeated 
or  to  change  the  element  parameters.  It  is  capable  of  carrying  out  dc 
analysis,  ac  analysis,  or  transient  analysis  of  a  circuit.  Furthermore, 
the  results  of  the  dc  analysis  can  be  directly  used  as  the  initial 
conditions  of  ac  or  transient  analysis  to  switch  into  ac  or  transient 
analysis. 
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I.  BASIC  CHARACTERISTICS 


Because  ECAL  uses  the  FORTRAN  language  to  compile  the  explanatory 
program,  it  can  be  operated  on  any  computer  with  a  Fortran  compiler. 

Its  symbol  group  is  a  subgroup  of  the  FORTRAN  symbols.  Its  data  for¬ 
mat  also  agrees  with  that  of  the  FORTRAN  language . 

The  symbols  are  formed  by  the  25  capital  alphabets  of  the  English 
language  and  the  10  Arabic  numerals  from  0  to  9  and  special  symbols, 
such  as  -+\  **— *  and  a  blank  space. 

The  statements  of  the  ECAL  language  are  also  composed  by  using 
cards  just  as  in  FORTRAN.  However,  each  statement  has  only  one  line 
which  corresponds  to  one  card.  There  is  no  continuation  card. 

The  foundation  of  electronic  circuit  analysis  of  the  ECAL  language 
is  the  nodal  point  method.  Between  nodal  points,  there  are  branches 
to  create  the  connections.  Each  branch  is  composed  of  a  sourceless 
element  and  several  independent  voltage  sources,  independent  current 
sources  and  correlated  current  sources.  Figure  1  shows  a  standard 
branch  circuit  formed  by  a  sourceless  element,  independent  power 
source  and  correlated  power  source. 

The  entire  analysis  is  carried  out 
on  the  basis  of  the  standard  branch. 

The  allowable  elements  are  resistors, 
capacitors,  inductors,  independent 
voltage  sources,  mutual  inductors, 
voltage  controlled  current  sources 
(abbreviated  as  mutual  conductance), 
current  controlled  current  sources 
(or  called  the  current  amplification 
coefficient),  semiconductor  diodes, 
semiconductor  triodes  and  nonlinear 
resistors. 

The  ECAL  language  uses  a  bending  line  method  to  describe  nonlinear 
elements.  The  nonlinear  resistor  N  can  be  used  to  describe  a  non- 


Fig.  1  Standard  branch 


linear  element,  such  as  a  tunnel  diode,  field  effect  transistor, 
voltage  stabilizing  transistor,  operational  amplifier,  etc. 

There  are  seven  output  formats  in  the  ECAL  language.  One  is  the 
standard  format.  When  the  format  is  not  specified,  this  format  is 
automatically  used.  It  tabulates  voltage  at  each  nodal  point,  voltage 
and  current  in  each  branch,  and  the  voltage  and  current  of  each  ele¬ 
ment  in  a  table  form.  If  it  is  an  ac  analysis,  it  also  provides  the 
power  of  each  branch  and  element.  Obviously,  this  method  may  involve  a 
huge  amount  of  output  data  especially  for  circuits  with  more  branches 
and  nodal  points  or  under  the  condition  that  the  analysis  frequency 
or  the  number  time  interval  points  is  high.  For  this  reason,  there 
are  three  more  table  output  formats  with  six  nodal  point  voltages, 
branch  current,  or  element  voltage  of  interest  as  well  as  three  curve 
output  formats  with  three  nodal  point  point,  branch  current,  or  ele¬ 
ment  voltage  of  interest. 

II.  APPLICATION  RANGE 

The  first  application  of  ECAL  is  dc  analysis.  It  can  be  used  to 
analyze  the  working  condition  of  the  various  serial  and  parallel  vol¬ 
tage  regulators  at  the  various  stages  of  different  dc ,  ac  amplifiers. 

It  can  also  be  used  to  analyze  the  amplification  factor  of  a  dc  ampli¬ 
fier  as  well  as  the  steady  state  voltage,  current  gains  of  the  medium 
frequency  equivalent  circuits  of  ac  or  pulsed  amplifiers. 

In  the  dc  analysis,  originally  it  does  not  treat  the  energy  stor¬ 
age  elements — capacitors  and  indicators.  However,  in  order  not  to 
destroy  the  completeness  of  the  circuit,  as  well  as  to  obtain  the  work¬ 
ing  points  from  dc  analysis  before  ac  or  transient  analysis,  therefore, 
inductors  and  capacitors  are  allowed  in  the  input.  In  the  computation, 
the  inductor  is  replaced  by  an  0.1  8  resistor,  while  the  capacitor  is 
replaced  by  a  10  M  ft  resistor  instead.  Consequently,  It  allows  the 
Input  of  the  dc  insulating  capacitors,  filtering  capacitors  or  compen¬ 
sating  inductors  of  a  multi-stage  amplifier  into  the  computer  without 
modification.  In  the  computation  of  dc  working  points,  the  explanatory 
program  automatically  will  use  the  aforementioned  resistance  values 

130 


instead.  The  computation  of  the  working  point  will  not  be  affected. 

Another  application  of  the  ECAL  language  is  ac  analysis.  It  is 
capable  of  conducting  frequency  characteristics  analysis  on  various 
wide  band,  narrow  band  amplifiers  and  filters  with  and  without  sources. 
For  the  AC  analysis,  the  explanatory  program  carries  out  the  analysis 
based  on  small  ac  signals.  It  is  usually  believed  that  nonlinear  prob¬ 
lems  do  not  exist  for  small  ac  signals  near  the  working  points. 
Therefore,  for  a  working  point  determined  in  the  dc  analysis,  if  the 
circuit  is  not  redescribed  after  switching  to  ac  analysis,  a  linear 
analysis  will  be  carried  out  for  the  corresponding  parameters  of  the 
nonlinear  elements  determined  by  that  working  point. 

One  important  statement  in  ac  analysis  is  the  frequency  range 
statement.  It  indicates  the  frequencies  to  be  analyzed.  Starting 
from  the  minimum  value,  it  is  increased  algebrically  or  geometrically. 
After  reaching  or  exceeding  the  maximum  value,  it  is  stopped.  Usually, 
for  a  narrower  analysis  range,  it  is  possible  to  use  an  equal  increment. 
Thus,  the  results  plotted  form  a  linear  coordinate.  If  a  wide  fre¬ 
quency  range  must  be  analyzed,  such  as  the  frequency  characteristics 
of  a  wide  band  amplifier  which  Is  usually  from  several  tens  Kz  to 
several  tens  mega  Hz,  then  it  is  not  suitable  to  use  a  linear  coordi¬ 
nate  or  equal  Increment  method.  At  this  time,  it  is  more  suitable  to 
use  the  common  ratio  increment  method.  Its  pattern  corresponds  to  a 
semi-log  coordinate. 

The  ECAL  language  can  also  be  used  to  conduct  transient  analysis 
on  circuits  and  systems.  Usually,  a  transient  analysis  Is  always 
done  with  respect  to  a  specified  input  signal  waveform.  For  example, 
the  output  waveform,  when  the  input  signal  is  an  impulse,  is  called 
the  impulse  response.  When  the  input  signal  is  a  unit  step  jump,  the 
output  is  called  a  unit  step  response.  These  are  very  Important  output 
results  in  transient  analyses.  In  addition,  sometimes  we  are  Interest¬ 
ed  in  the  response  of  the  system  or  circuit  to  a  certain  input  waveform. 
For  example,  a  ramp  waveform  (linearly  increasing  waveform)  and  expon¬ 
entially  increasing  waveform,  etc.,  are  commonly  used.  The  ECAL  lan¬ 
guage  provides  the  feasibility  to  describe  the  input  signal  which  is 
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also  comparatively  simple.  It  involves  the  use  of  an  input  state¬ 
ment  INPT.  It  uses  the  positions  of  four  coordinate  points  to  provide 
the  variation  of  the  input  signal.  Thus,  it  is  possible  to  describe 
a  step,  a  ramp,  a  square  wave,  a  trianguar  wave,  a  sawtooth  wave  and 
a  trapezoidal  wave.  It  is  also  possible  to  use  three  segments  of 
straight  lines  to  replace  an  exponential  or  logarithmic  curve  approx¬ 
imately. 

III.  PROGRAM  PLOW  CHART 

In  order  to  suit  the  batch  process  method  [5]  of  the  operating 
system  of  the  Felix  C-256  computer  presently  in  use,  the  explanatory 
program  of  the  ECAL  language  has  also  adopted  a  batch  process  form, 
i.e.,  several  circuit  analysis  jobs  can  be  processed  at  the  same  time. 
Each  job  is  not  mutually  correlated.  A  job  not  only  can  perform  a 
dc  analysis,  but  also  an  ac  or  transient  analysis.  It  is  also  possible 
to  change  the  element  parameters  repeatedly.  This  format  can  be  con¬ 
verted  into  the  man-machine  dialog  format  easily  with  user  teletype 
terminals  or  CRT-keyboard  terminals. 

The  explanatory  program  can  be  approximately  divided  into  three 
parts:  the  main  program — reading  in  the  headings  and  the  electronic 

circuit  descriptive  statements,  filling  out  the  table  of  elements  and 
controlling  the  process  of  analysis;  and  two  subroutines — calculating 
and  providing  output  results  of  the  dc,  ac  and  transient  analyses  of 
the  circuit,  respectively. 

The  program  flow  chart  is  shown  in  Figure  2. 

IV.  EXAMPLES  OF  THE  ANALYSIS 

A  low  frequency  amplifier,  as  shown  in  Figure  3»  has  22  branch 
numbers  and  nine  nodal  point  numbers.  The  nodal  point  number  is  circled 
to  distinguish  from  the  branch  number.  The  18th  branch  resistance  is  a 
negative  feedback  resistance.  The  circuit  analysis  program  is  shown 
In  Table  1.  The  first  execution  statement  EXEC  carries  out  the  dc 
analysis  to  calculate  the  working  point.  Then,  the  dc  power  supply 


KEY  TO  FIGURE  2  (page  133): 

1 — start;  2 — initial  job  statistics;  3 — read  headings;  4 — is  it  a 
blank  card;  5 — yes;  6 — no;  7 — yes;  8 — is  the  first  part  of  the  special 
identification  code  ENDJ ;  9 — no;  10 — PRINT  NJOB  and  headings;  11 — 
clear  working  area  and  element  table;  12 — read  one  card;  13 — is  it 
ENDJ;  14 — yes;  15 — stop;  16 — yes;  17 — is  it  FNSH;  18 — no;  19 — is  it 
TIME;  20— no;  21  — is  it  FREQ;  22— no;  23— is  it  PRNT;  24— no;  25— 

Is  it  PLOT;  26 — no;  27 — is  it  INPT;  28 — no;  29 — yes;  30 — yes;  31 — yes; 
32 — yes;  33 — yes;  34 — note  time  parameter;  35 — note  frequency  para¬ 
meter;  36 — note  print  type;  37 — note  plot  type;  38 — note  input  signal 
oarameter;  39 — is  it  MODF;  40 — no;  41 — is  it  four  blank  spaces; 

42 — no;  43 — is  it  EXEC;  44 — no;  45 — none  of  the  above.  This  card  is 
wrong.  PRINT  "THIS  CARD  IS  WRONG";  46 — yes;  47 — yes;  48 — yes; 

49 — note  number  of  times  of  value  charges;  50 — enter  the  elements 
according  to  the  type  into  the  element  table;  51 — check  for  mistakes 
in  the  input  statements  and  the  correlation  matrix;  52 — yes;  53 — 
situation  when  printing  error  occurs;  54 — fill  the  matrice  A;  55 — no; 
56 — yes;  57 — no;  58 — no;  59 — yes;  60 — call  subroutine  DTANR ;  61 — 
call  subroutine  ACANR 


is  removed  and  a  5mV  ac  voltage  signal  is  added.  At  R-^g  =  10  MD , 
it  corresponds  to  the  situation  that  the  negative  feedback  linkage  is 
broken  off  to  obtain  the  frequency  characteristic  curves.  Then,  R^g 
=  20  Kft  is  used  to  obtain  the  frequency  characteristic  curve  with  the 
feedback.  The  range  of  analysis  is  from  1  Hz  to  1  mega  Hz. 


The  circuit  for  a  single-shot  trigger  is  shown  in  Figure  4.  The 

t  h 

trigger  pulse  Is  added  to  the  21  branch.  Different  triggering  capa¬ 
citance  has  an  effect  on  the  output  waveform.  Here,  a  value  change 
statement  is  used  to  change  The  response  curves  corresponding  to 

C13  =  4000PF,  8000  PF  and  12,000  PF  were  calculated.  The  trigger 
pulse  is  a  square  wave,  20  microseconds  wide.  The  amplitude  is  four 
volts.  The  program  is  listed  In  Table  2. 
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ELECTRONIC  CIRCUIT  ANALYSIS  LANGUAGE  (ECAL) 

Chert  Chenghang 

(Northwestern  Polyicchnical  University) 

Abstract 

The  computer-aided  design  (CAD)  technique  is  one  of  important  develop¬ 
ments  in  computer  application.  CAD  technique  has  already  become  an  important 
branch  of  computer  science.  Special  language  for  electronic  circuit  analysis  is  the 
foundation  of  computer-aided  circuit  design  (CACD)  and/or  computer-aided  cir¬ 
cuit  analysis  (CACA). 

Electronic  circuit  analysis  language  (ECAL),  a  special  language  for  circuit 
analysis*  is  comparatively  simple  and  convenient  for  engineering  application. 
Statements  of  ECAL  are  executed  explanatorily  by  FORTRAN  language.  So  far, 
ECAL  can  be  used  to  make  DC*  AC  and  transient  analysis  of  a  circuit*  and 
the  results  of  DC  analysis  can  be  regarded  directly  as  initial  conditions  for  AC 
and  transient  analysis. 

Both  linear  and  nonlinear  elements  can  be  taken  into  account  in  ECAL. 
Linear  elements  may  consist  of  resistors,  capacitors,  inductors,  mutual  inductors* 
independent  current  sources,  independent  voltage  sources  and  dependent  current 
sources  under  voltage  control  or  current  control.  Nonlinear  elements  m.-y  include 
diodes*  transistors  and  nonlinear  resistors.  Therefore,  this  language  is  suitable 
to  the  analysis  of  linear  circuits  or  systems  as  well  as  nonlinear  ones. 


AN  EXPERIMENTAL  INTERACTIVE  COMPUTER  GRAPHICS  SYSTEMS 
FOR  FREE-FORM  SURFACE  DESIGN 


Zheng  Hulling,  Wang  Zhisheng,  Lu  Hongjia  and  He  Tianbao* 
Shanghai  Aircraft  Manufacturing  Factory 


ABSTRACT 

An  interactive  graphics  system  for  free-form  surface 
design  is  introduced  in  this  paper.  The  system  was  esta¬ 
blished  on  the  basis  of  the  cubic  uniform  B-spline  theory. 

In  order  to  obtain  better  simulated  results,  the  basic  func¬ 
tion  with  a  quadruple  knot  at  both  ends  was  chosen.  Its 
major  functions  are  as  follows: 

1.  It  displays  the  three-dimensional  model  on  a  screen 
or  a  plotter.  It  also  allows  the  three-dimensional  coordi¬ 
nate  system  to  undergo  real  time  transformation  with  respect 
to  the  model. 

2.  Through  the  use  of  a  method  involving  the  fairing 

of  the  curvatures  of  discrete  points  sequentially,  the  exter¬ 
nal  shape  of  the  model  is  faired. 

3.  It  provides  two  local  modification  methods  for  the 
model  surface  and  also  ensures  C2  continuity. 

4.  It  permits  the  display  of  any  arbitrary  cross-section 
of  the  model. 

A  computer  aided  shape  design  interactive  system  was 
established  using  a  DJS-6  computer  and  a  model  "751"  optical 
pen  graphics  display  which  was  developed  and  constructed  by 
Sian  JIaotung  University.  Presently,  it  is  taking  its  shape. 
However,  perfection  and  further  development  are  continuing. 

The  B-splines  method  was  chosen  to  form  the  curve  (surface) 
because  of  considerations  in  geometric  intuition,  ease  of  com¬ 
putation  and  ease  of  control  and  modification.  It  is  graphi¬ 
cally  called  the  "characteristic  ployhedron  surface  design 
method".  A  brief  introduction  to  the  major  functions,  con¬ 
clusions  and  application  examples  was  presented  in  this  paper. 
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I.  THE  STRUCTURE  AND  DISPLAY  OF  THE  CURVE  (SURFACE) 

With  regard  to  the  design  of  external  shape,  the  use  of  a  cubic 
uniform  B-splines  method  is  very  suitable.  In  order  to  improve  the 
end  point  characteristics,  a  basic  function  [1]  with  a  quadruple  knot 
at  each  end  was  used  (Figure  1).  In  comparison  to  a  usual  cubic 
uniform  B-splines  curve,  a  difference  only  occurs  in  the  two  segment 
of  curves  on  the  ends.  In  the  middle,  they  coincided  completely 
(Figure  2).  It  has  very  good  extrapolation  conditions  on  the  front 
end  which  satisfy 


Fif .  1  B-apHnri  with  a  quadruple  kaot  at  each  tad  poiat 


The  rear  end  is  similar  to  the  front  end.  From  equation  (1)  we  can 
see  that  the  curve  passes  through  the  front  end  point.  Furthermore, 
it  is  tangent  to  the  front  (real)  end.  From  equation  (2)  we  can  see 
that  usually  0 )  is  not  equal  to  zero.  Therefore,  the  method  of 
using  polygon  equi-distance  extension  to  ensure  that  the  usual  cubic 
uniform  B-spline  passing  through  the  end  points  is  superior. 


Fi(.  2  Tba  cormpoadinl  B-»plin*  eurra 

Key:  1 — vortex;  2 — knot;  3 — the  usual  uniform  B-spline 
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Prom  equation  (1)  we  can  also  see  that  the  conditions  to  be 
satisfied  at  the  end  point  are  consistent  with  those  of  the  Bezier 
curve.  Especially  when  only  one  segment  of  curve  exists  (four  vert¬ 
ices),  this  type  of  curve  rigorously  deteriorated  into  the  Bezier 
curve.  Therefore,  this  system  permits  the  combined  use  of  both  struc 
turing  methods. 

For  the  convenience  in  use,  a  matrix  expression  of  this  type  of 
curves  under  various  conditions  is  derived  from  the  deBoor-Cox  intera 
tion  equation.  For  example,  when  «>8,  the  first  two  segments  of 
curves  can  be  expressed  as: 
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The  last  two  expressions  are 
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The  middle  segments  are  consistent  with  the  usual  cubic  uniform  con¬ 
dition 
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Summarizing  the  coefficient  matrices  under  various  possible  con¬ 
ditions  (not  limited  by  m  and  i),  we  discovered  that  there  are  only 
18  different  columns.  In  the  program,  these  18  columns  are  placed  in 
a  4xl8  group.  Each  time  four  columns  are  extracted  based  on  the 
various  m  and  i  values  to  form  the  coefficient  matrix  of  the  B-spline 
of  this  segment.  Thus,  the  number  of  operation  is  reduced  and  the 
internal  storage  capacity  is  saved.  This  4x18  numeral  group  is 
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As  for  the  curved  surface,  the  range  of  boundary  multiple  knot 
effect  is  shown  as  the  shaded  area  in  Figure  3.  So  long  as  the  curve 
is  treated  well,  the  curved  surface  problem  is  also  naturally  resolved. 

Figure  4  is  the  characteristic  polyhedron  of  an  aircraft  dis¬ 
played  by  using  the  VVDP  light  button  through  the  optical  pt.i  points. 
The  point  curved  surface  display  light  button  SUDP  can  display  an 
arbitrary  projection  of  the  curved  surface. 
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Fi|.  3  The  areas  influenced  by  the  multiple  knots  on  a  surface 


ffl  *  03  5  TLSlftiiSM** 

Fi*.  4  The  characteristic  polyhedron  of  an  aircraft  Fi*.  5  Display  of  the  surface  of  an  aircraft 


II.  SOLVING  FOR  THE  VERTEX  AND  DISPLAYING  THE 

CHARACTERISTIC  POLYHEDRON 

In  computer  aided  design  and  manufacturing,  the  composite  prob¬ 
lem  is  frequently  encountered.  It  is  necessary  to  start  from  the 
known,  discrete  vlaue  point  to  solve  for  the  vertex  of  the  character¬ 
istic  polyhedron.  Subsequently,  the  B-spline  curve  passing  through 
the  above  mentioned  point  Is  created. 

This  system  considered  the  vertex  problem  in  depth,  Including 
various  shape  value  point  numbers  n(the  effect  of  end  point  is  diff¬ 
erent),  the  three  end  points  condition — free  end,  the  usual  cubic  uni 
form  B-splines  of  the  known  tangential  vector  or  equal  distance  exten 
sion  and  various  combinations  of  them.  There  are  26  situations  in 
total.  For  exmaple,  when  n  >  5,  the  computation  formula  in  solving 
for  the  vertex  under  the  condition  that  both  ends  are  free  ends  is 


1H2 
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In  equation  (7),  the  second  row  and  the  second  last  row  come 
from  the  knot  condition  (choosing  *  =  3(3))  of  the  first  and  last 
segments.  At  this  time,  m  *  n.  Let  us  take  the  tan¬ 

gential  vectors  and  of>  both  ends,  then  the  two  rows  should  be 
changed  to: 

—  3v,+3k,  =  5„  -3Vm.t  +  Zv„~Ut  (Q) 

In  addition,  when  m  *  n  +  2,  which  corresponds  to  the  equal  dis¬ 
tance  extension  condition,  the  first  two  rows  should  be  changed  to: 

vWi.  V 1  ~ 2 V:  +  V .1  -  0  ^ 

The  last  two  rows  are  similar  with  m  =  n  +  2. 

The  vertex  of  a  curved  surface  can  be  obtained  through  the  shape 
value  point  after  two  such  processes. 

92 

This  system  can  also  assure  that  the  B-spline  obtained  from  the 

vertex  contains  a  pre-determined  straight  line  segment,  including  the 

2 

pre-assigned  starting  and  end  points.  Furthermore,  a  C  continuity 
o 

exists  there  (C  is  a  second  order  continuity).  For  a  curved  surface, 
this  means  that  it  is  allowed  to  contain  a  pre-appointed  plane. 
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III.  PARTIAL  MODIFICATION  OF  THE  CURVE  (SURFACE) 

AND  THE  DISPLAY  OF  THE  MODIFIED  POLYHEDRON 

This  is  a  power  means  of  a  computer  aided  external  shape  design 
system  with  man-to-machine  dialog  capability.  It  should  be  compara¬ 
tively  more  effective  in  solving  this  problem  by  using  the  partial 
support  characteristics  of  the  B-splines. 

Two  types  of  partial  modification  capabilities  were  arranged  in 
this  system: 

1.  Direct  modification  of  the  vertex  of  the  polyhedron.  Through 
the  optical  pen  and  the  keyboard  to  send  in  the  corresponding  modifi¬ 
cation  information — we  can  change  the  vertex  number,  point  number  and 
its  coordinate  value  until  the  shape  is  considered  to  be  satisfactory. 

2.  Modification  of  the  vertices  of  the  polyhedron  by  using  a  new 
curve  (surface)  value.  It  is  required  that  the  new  curve  surface 

generated  after  the  modification  of  the  vertices  must  pass  through  the 

2 

newly  appointed  points.  Furthermore,  C  continuity  must  still  be 
ensured. 

Because  basically  cubic  uniform  B-splines  are  used,  the  principle 
of  modification  is  very  simple.  Using  the  curve  modification  in  Figure 
8  as  an  example,  ?cf'  are  the  new  given  points.  Our  treatment 

is  to  discard  the  shape  value  control  at  B  and  F  (which  is  frequently 
near  a  round  angle  transition;  therefore,  it  is  acceptable).  The  new 
vertices  Vc*%  vD\  v»f,  are  solved  simultaneously. 

2  —  j 

~y~Vc'+  ~g-V»'  =  ?c, - i-p, 

~g-  Vc'  +  -y- vV  +  y-vE/  =  p0,  (9) 

~~Vo>  +  j -JV-p,, — Lp, 
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Fig.  8  Determination  of  the  *er  tiers  bj  new  giren  points 
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Because  only  three  vertices  Vc',Vp'.vb'  are  changed,  therefore, 

the  new  curve  determined  by  the  polyhedron  . VsViVc'Vd'Vi'VfVg . 

2  _  _ 

can  satisfy  C  continuity.  Only  in  the  P*Pc  segment,  it  deviates 
from  the  original  curve.  Furthermore,  it  rigorously  passes  through 
Pc\Pd\Pb'  ;  therefore,  it  agrees  with  the  modification  requirements. 
Moreover,  the  computational  work  load  is  very  small.  Similarly,  it 
can  be  used  on  a  curved  surface.  The  corresponding  surface  knots  can 
also  be  modified  locally  in  order  to  save  the  computational  time.  It 
should  be  pointed  out  that  the  disadvantage  of  this  treatment  is  that 
the  convex  boundary  shape  under  surface  modification  cannot  be  satis¬ 
fied  (the  modification  zone  is  limited  to  above  the  rectangular  region 
of  the  parametric  plane).  We  computed  an  example  which  involved  the 
insertion  of  a  large  hemisphere  on  a  wing.  After  modifying  20  points, 
including  the  transition  region  between  the  hemisphere  and  the  wing, 
we  obtained  very  good  results  (Figure  9). 
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IV .  CURVE  FAIRING 


The  fairing  of  the  B-splines  can  be  reduced  to  the  fairing  of 
the  corresponding  characteristic  polyhedron  vertex  sequence  [2].  We 
adopted  the  curvature  method  introduced  in  that  book  to  carry  out 
fairing  with  specific  reference  to  the  bad  points  of  the  vertex  seq¬ 
uence  of  the  polyhedron.  Its  basic  principle  is  that  the  ordinate  of 
the  vertex  must  be  modified  in  order  to  ensure  the  quadruple  differ¬ 
ence  of  the  ordinate  with  respect  to  the  arc  length  at  that  vertex  is 
zero  (A 0  -)» 


The  so-called  curvature  of  point  p.,  is 

circle  passing  through  the  three  points  p,M,  p,„ 

AT  —  ‘•A/  _  _  2sin<P, 

'  Uuji  ~  i:  " 


the  curvature  of  the 

(10) 


Fif.  10  Tbc  curvature  of  >  tpecifie  circle  of  Pi 


where  A.  is  the  algebraic  area  of  the  triangle  formed  by  p  .,p.p,,, 
Assuming  that  P4  is  determined  to  be  a  bad  point  in  Figure  11,  then 
the  modification  quantity  has  the  following  computational  formula: 


Equation  (11)  was  described  in  [2].  Its  detailed  proof,  program 
flow  chart  and  effectiveness  were  given  in  [ 5 ] -  Prom  the  effective¬ 
ness,  it  is  comparatively  ideal.  Furthermore,  because  the  computation 
al  process  does  not  Involve  the  superpositioning  of  curves,  the  fair¬ 
ing  speed  is  fast.  As  for  a  spatial  curve,  it  is  converted  into  two 
projection  curves  to  be  treated  separately. 

V.  THE  COMPUTATION  AND  SECTION  DISPLAY  OP  AN  ARBITRARY 
CROSS-SECTION 

In  order  to  rigorously  control  the  external  shape,  especially 
some  key  positions,  this  function  is  also  mandatory.  Due  to  the 
inclination  of  the  sectional  coordinate  system  with  respect  to  the  sur 
face  coordinate  system,  many  complicated  details  are  brought  into  the 
computation  of  the  section.  In  this  system,  the  vertices  and  knots 
of  the  polyhedron  are  transformed  into  the  sectional  coordinate  system 
based  on  the  geometric  invariance  of  the  B-spline  surface.  Subse¬ 
quently,  the  computation  of  the  inclined  sectional  external  shape  is 
transformed  into  the  normal  section  (5-0)  .  Therefore,  it  is  much 

faster.  Furthermore,  it  is  more  convenient  for  the  digital  con¬ 
trolled  plotting  of  the  sectional  shape,  computation  of  the  incline 
angle  and  digitally  controlled  shape  modifications.  After  the  compu¬ 
tation,  the  system  can  automatically  transform  the  vertices  and  knots 
back  to  the  surface  coordinate  system  to  be  used  for  the  next  sec¬ 
tional  computation.  The  sectional  plane  parameters  are  sent  into  the 
system  by  the  keyboard.  The  step  length  of  the  extrapolation  on  the 
section  line  can  provide  the  A*  of  the  sectional  coordinate  system  or 
the  parametric  Increments  AU  or  AW.  It  is  selected  by  the  user. 

Figure  12  shows  the  inclined  section  of  a  wing  bulge  (see  Figure 
9).  The  step  length  chosen  is  AW  =  0.2. 

VI .  IDENTIFIED  PROBLEMS 

In  addition  to  the  aforementioned  light  buttons,  functional  keys, 
the  present  system  also  includes  over  20  optical  buttons  and  function 
keys  such  as  picture  change  (magnification,  reduction,  displacement, 


9 


rai2 

Fi*  12  The  section  in  *  sectional  coordinate  svstrni 

rotation),  automatic  arrangement  of  the  picture,  graphics  output  and 
three-dimensional  coordinate  system  transformation  and  three-dimen¬ 
sional  surface  symmetry,  etc.  The  function  mentioned  above  cannot 
meet  che  actual  requirements  yet.  There  is  a  lot  of  work  to  be  done. 
For  example,  they  include  the  further  modification  of  the  GSP  software 
system,  the  establishment  of  a  fast  algorithm  for  the  display  of  curve 
(surface),  the  formation  of  a  transition  surface,  the  continuity  prob¬ 
lem  of  the  curve  (surface),  the  surface  plate  problem  of  a  non-quadri¬ 
lateral  surface ...  and  other  applied  programs,  such  as  the  connection 
between  the  strength  and  aerodynamic  computation  programs. 

In  addition,  the  present  hardware  conditions  are  far  less  than 
those  required  to  satisfy  the  requirements. 
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AN  EXPERIMENTAL  INTERACTIVE  COMPUTER  GRAPHICS 
SYSTEM  FOR  FREE-FORM  SURFACE  DESIGN 

Zheng  Huiling,  Wang  Zhisheng,  Lu  Hongjia,  He  Tianbao 

C Shanghai  Aircraft  Manufacturing  Factory ) 

Abstract 

An  interactive  computer  graphics  system  for  free-form  surface  design  is 
described  in  this  paper.  and  application  of  cubic  uniform  B-splines  is  suggested 
as  the  fundamental  method  of  surface  modeling.  In  order  to  obtain  a  better  appr¬ 
oximation,  the  basic  function  with  a  quadruple  knot  at  each  end  point  is  utilized. 

In  respect  of  curve  construction,  there  are  only  18  columns  in  each  of  va¬ 
rious  coefficient  matrices  for  the  different  numbers  of  vertexes  and  the  different 
sequences  of  curve  segments,  and  the  matrix  expressions  of  curves  are  also  gi¬ 
ven.  As  a  result,  real  benefit  is  gained  for  reducing  the  storage  capacity  and 
increasing  the  computational  speed.  In  the  inverse  calculation  for  the  vertexes 
26  various  composite  cases  are  summarized,  which  permit  to  maintain  the  ori¬ 
ginal  straight  line  segments. 

It  is  proposed  to  fair  the  given  curves  by  fairing  the  curvature  sequence 
of  the  discrete  vertexes,  and  preliminary  practical  experience  is  also  given. 

This  paper  offers  a  simple  and  fast  engineering  algorithm  for  modifying  of 
the  surface  data  in  a  local  region,  while  the  surface  still  remains  in  second 
order  continuity. 

In  the  experimental  system  there  are  20  functional  buttons  (embracing  the 
functions  of  coordinate  transform.  3— dimensional  symmetry,  construction  of  an 
arbitrary  section  etc.  ).  A  success  has  been  made  in  surface  modeling  of  aircraft, 
automobile  and  ship. 
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BEZIER  PLOTTING  THEOREM  AND  GEOMETRIC  CHARACTERISTICS 
OF  CUBIC  BEZIER  CURVES 


Shi  Fazhong  and  Wu  Junbeng  * 

Beijing  Institute  of  Aeronautics  and  Astronautics 


ABSTRACT 

The  geometric  characteristics  of  the  Bezier  curve 
have  been  studied  in  depth  by  P.  E.  Bezier  using  the  fast 
end  curve  [1]  and  by  Su  Boqing  and  Liu  Dengyuan  using 
simulated  projection  transformation  [3-5]-  The  method 
presented  in  [1]  by  Bezier  to  find  the  points  and  their 
tangents  on  the  Bezier  curve  using  geometric  plotting  is 
the  Bezier  plotting  theorem.  This  paper  analyzed  the  geo 
metric  characteristics  of  plane  cubic  Bezier  curves  based 
on  the  plotting  theorem.  It  pointed  out  that  the  X,  y 
(or  I,  y)  shown  in  Figure  2  are  a  pair  of  invariant  quan¬ 
tities  determining  the  geometric  characteristics.  The 
complete  plane  diagram  of  X,  y(or  X,  y)  was  given  (see 
Figure  3).  Some  geometric  characteristics  of  the  spatial 
cubic  Bezier  curves  were  discussed. 


I.  THE  DERIVATION  OF  THE  BEZIER  PLOTTING  THEOREM 

f"  Vi 

The  Bernstein  expression  of  the  n  order  Bezier  curve  is 

* 

“)«  2  «•*(■)•?/  ««co,  n 

i  »  o 

$..,(«  )  =  c'«'(  1  -  «)•*' 

where  is  the  cusp  vector  of  the  Bezier  characteristic  polygon,  9.„(  « ) 

fy  h 

is  the  n  Bernstein  basic  function. 

The  geometric  plotting  method  presented  by  Bezier  in  [1]  to  find 
the  points  and  their  tangent  on  the  curve  is  shown  in  Figure  1.  We 
called  it  the  Bezier  plotting  theorem.  It  can  be  expressed  in  terms  of 
the  following  set  of  iteration  equations: 

»  -  p  (2a ) 

P(u)^  2  "  )•$»' ’’(  "  )  A  =  o.i,  •••,«! 

i  * » 


(la) 

(lb) 
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S/'  («0  =  s7'-,;(«0+  u  (  U)  (2b) 

S, c*'(«  )«S,  f  2c) 

p'  (  u  )-  n  S.c-,5(«)S7r,3(  «  )  <2d) 

p*(  u  )=  n  (  it  -  i  )CS1c-‘3(u)57"i)  (  «  )  -S0c-*:(u)S,c-‘;(  u  »  (2e) 

Chang  Gengzhe  and  Wu  Junhen  first  provided  the  proof  [2].  The  follow¬ 
ing  can  also  be  derived  directly. 


Fi*  1  Application  of  the  plotting  theorem  to  determining  point* 
and  their  tangent*  on  the  Ber.ier  curve 


Using  ei’*ci.,+ci:1,  and  changing  the  index  j,  we  can  rewrite  (la) 
into  "-1 


?( « )-  2  ci.,u'(  1  -  u  )‘-'s,+  2  c':{  «'( I  -  U  )•-'£, 

/  »  0  j  rn  l 

"  -  l  a  -  i 

■  2  c>-iu>(  1  ~  "  1  -  “  )S;  +  «S/*i)“  2  «*-..»(  “  )S/‘3(  “  ) 


where  «)-(!-■  )s,+«Sm- s,-+  «  Cs,.,-Si) 


Repeating  the  above  process,  we  know  that  (2a)-(2c)  are  valid. 
Using  (lb)  to  find  its  partial  derivatives  with  respect  to  u,  we  get 

</:,(.  »  )  -  i »  C*.t.,.,(  »  »  )) 

g':>(  «)-»(»-  1  )Cp«.„/-»(  “  »  )  +  *.»./(  “ 

Again-,  finding  the  derivatives  of  (la)  with  respect  to  u,  and  using 
the  two  above  equations,  we  know  (2d)  and  (2e)  are  valid.  The  higher 
order  derivative  vectors  can  be  derived  by  further  steps. 
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II.  GEOMETRIC  CHARACTERISTICS  OF  PLAN/ R  CUBIC 
BEZIER  CURVES 


The  eauation  of  a  plane  cubic  Bezier  curve  is 

P(  “  )«(  1  -  “  )*S,  +  3u  (  1  -  u  ySi  +  Zu'C  1  -  u  )s,  +  u*S,  u«(0,l3  (3) 

where  S/.  /  *  0,  1 ,  2,  3  %  «/■* si—  S/m“  (*/.  y>»  0 ),  j  >■  1 ,  2,  3  correspond  to  the 
four  co-planar  vertices  and  three  vectors  of  the  characteristic  tri¬ 
angle,  respectively.  If  the  lines  coincide,  then  the  curve  becomes 
a  straight  line.  No  discussion  will  be  necessary. 

The  geometric  characteristics  of  a  plane  cubic  Bezier  curve 
include  whether  a  single  point  (a  cusp  or  a  double  point)  or  an  in¬ 
flexion  point  (an  inflexion  point  or  two  inflexion  points)  exists  or 
whether  the  curve  is  convex  or  not. 

We  used  the  two  values  X,y  as  shown  in  Figure  2  to  express  the 
geometric  characteristics  of  a  plane  cubic  Bezier  curve.  Once  the 
observed  X,y  are  determined,  the  geometric  characteristics  of  the  curve 
are  determined.  It  Is  not  related  to  the  amplitude  and  direction  of 
the  vector  of  the  sides.  This  Indicates  that  X,u  are  a  pair  of  invar¬ 
iant  quantities  or  geometric  characteristic  control  parameters  which 
determine  the  geometric  characteristics  of  the  plane  cubic  Bezier 


curve , 


C o/,  5^3 1 


«■  C/S,,  0.3,  then 


we  get 


y>  y> 


i  M  d _  A  A  A  _  _  A _ 

A,t  1  -  x  ’  ”  y  ’  15  mT  -  x ) 


We  are  designating  the  parameters  corresponding  to  the  appearance 

of  an  Inflexion  point  and  a  cusp  on  the  plane  cubic  Bezier  curve  as 

and  u  ,  respectively.  The  two  parameters  corresponding  to  the  double 
c 

point  are  noted  as  u^_  and  u2»  Then,  the  equations  of  inflexion  point, 
cusp  and  double  point  of  the  curve  and  their  corresponding  equations 
of  the  plotting  theorem  are  the  following: 


iw'^V  — 


w&Sm 


r 


ra  2  fttEHj&UtoX,  I*  XfX,  ? 

Fig.  2  The  X  ,  M  and  X  ,  I*  of  a  characteristic  trilateral 


X 


+  1  , 


(5) 


inflexion  point  equation?/(n,)x  p'(ti,)«o,  p'(u, )?*o.m(  o,  i  )=»  Stm(u,)Stcn  («>) 

xS‘,,(«l)S1(l‘(«.)-5,  f,‘u(«r)**$itM(*)  (6) 

cusp  equation  p' (lie)- o,  m(  0 ,  1  )==^s.U3(««c)*St,3(«c)  (7) 

double  point  equation 

#(«,)- p(«h).  o<»I<ii1<i*o<iiI<iiI<i=»s.t*’(«»)-‘S.,'*3(“»)  (8) 

where  the  parameter  region  is  unchanged.  It  is  discussed  separately 
in  the  following: 


1. 


Inflexion  point  equation:  From  equation  (6),  we  can  get 
(Alt  +  Ati— Ats)u*  +  (yiu— ■  2Alt)ui  +  Altm  0 

Using  equation  (4),  the  above  equation  can  be  rewritten  as: 

(2  +  H-X)u?~(l  +24M+H-0  (9) 

whenA-P  +  2,  we  have  s,1,,(«,)$l;,J  («,)  -  Slt,,(n,)Slt,J  («,) 

The  curve  becomes  a  single  inflexion  segment  [6].  It  Is  possible  to 
conduct  a  simulated  projection  to  transform  It  into  a  usual  cubic 
polynomlnal  [3].  Equation  (9)  has  only  one  root: 

** 


1  +  2H 


(10) 


In  order  to  let  ut«(0,l)  ,  It  is  necessary  to  have 11  >o  ori*<-l 
The  curve  has  a  single  inflexion  point  at  (0.1)  as  shown  in  Figure 
example  1. 


When  \  +  »  +  2 


then 


,1+2^  ±  \/  4  n  (  x  -  i )  + i_ 

'  2(2  +  ■»->.) 


(11) 
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The  determinant  £=»4H(X  —  l )+ 1  has  the  following  several 
situations : 

(1)  A  >  0,  there  are  two  real  roots  v, ,  however,  they  may 
not  all  be  on  (0,1). 

1)  If  n>0.X>iorn<0.X<i  ,  then: 


(  X  >  2  +  H 

«/,< 

when 

a. 

]  1  <  X  <  2  +  M 

at 

«/;> 

'  X  -  1 

U,.  = 

when 

M  =  0 

.  X  <  1  It.f,  ui j 

-  o. 

(  x  <  2  +  H 

when 

H  <  fl 

• 

at 

2  +  H<X< 1 


«/;<  0  J 
“»i>  1  . 


1C 


i.e.,  there  is  always  a  root  outside  (0,1).  There  is  only  one  inflex¬ 
ion  point  on  (0,1)  as  shown  in  the  diagrammatic  examples  2-4.  As  a 
special  case,  when  a,//  «,  X ,  n  =oor£-c^,  as  shown  in  Example  5.  Only 
at  this  time  the  inflexion  point  parameters  cannot  be  expressed  by 
X,  u.  It  is  possible  to  assume  “j**  k  a„k  >  o,  then  we  can  obtain  the 
following  from  equation  (6): 


Ufs 


-  1  +y/  k 

k  -  "l 

1/2, 


h¥*  1 
h  =  1 


(12) 


We  define  that  after  putting  them  in  sequence,  the  vector  of  one 
side  has  a  rotational  angle  with  respect  to  the  previous  one  of  the 
same  sign  and  its  absolute  value  is  always  less  than  180° ;  the  plane 
polygon  has  an  unchanged  direction  of  rotation.  Otherwise,  its  direct¬ 
ion  of  rotation  is  changing.  Then  by  summarizing  the  situation  that 
the  curve  has  an  inflexion  point  on  (0,1),  we  find  the  necessary 
and  sufficient  conditions  for  a  plane  cubic  Bezier  curve  to  have  an 
inflexion  point  on  (0,1)  are  that :  v  >  0 ,  X >  ior  ^ <  o ,  X <  l ,  i.e.,  the 

corresponding  characteristic  trilateral  changes  its  direction  of  rota¬ 
tion. 


-  v  ~- 
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2)  If  n>o,  i j^j  ,  then  the  two  real  roots  ur, .  mt<(  o ,  l ), 

i.e.,  the  two  Inflexion  points  (also  called  the  double  inflexion 
point)  must  appear  on  (0,1).  As  a  special  case,  when  ScSt//S,S,. 

o<^»x<l  ,  then  1/2  .  When  the  extension  of 

the  sides  SqS2  and  intersects,  i.e.,  when  0O<X<i,  there  must 

be  a  double  inflexion  point  as  shown  in  example  6.  ^ >  o ,  l  -  4^-  <  X <  1 

are  the  necessary  and  sufficient  conditions  for  the  plane  cubic  Bezier 
curve  to  have  two  inflexion  points  on  (0,1). 


(2)  When  A  =  0,  there  is  a  double  root.  From  the  following  dis¬ 
cussion  on  the  cusp  equation,  we  know  that  this  is  the 

?'(«c)  =  0=4>  S.'-l\uc)=slilj(uc) 

condition.  The  two  inflexion  points  coincide  with  each  other  which  means 
it  is  a  cusp  as  shown  in  examples  7-8.  Therefore,  we  have 

X  -  1  -  A  #  **>o 

44  (13) 

^ y2ti ,  u>o 

The  necessary  and  sufficient  conditions  for  the  plane  cubic  Bezier 

curve  to  have  a  cusp  on  (0,1)  are:  w  1 

n  >  o ,  x  =  i  --jv* 

(3)  A<  0  x  <  i n>o  »  there  is  no  real  inflexion  point. 
From  the  following  discussion  on  the  double  point  equation,  we  know 
that  when  A  is  above  a  certain  lower' limit ,  the  curve  will  have  a 
double  point  on  L0,1)  or  (0,1]. 

2.  Cusp  equation.  Here  we  will  directly  find  the  cusp  from  the 
cusp  equation  (7).  From  equation  (7),  we  get: 

SiC,,0«cW”  (“c)  «  ~ ~  S,l,i(uc)Sf»  (uc) 

=»(o,-23t+3;!)u=+  2  (S.-a.X.  +  a.-o  ) 

We  can  obtain  X  *  l  «c*  - ^>o. which  is  consistent  with  (13)  • 

4  1  x  2  H 


3.  Double  point  equation.  From  equation  (8),  we  get 


we  can  obtain 


(a,~23,  +  3,)(i4j+!ilu.  +  ii*,)+  3  (3. -a, )(«,+«,)  +  33,-0 


t+2**  ±  »  12“  (  1  -  x  )- 3 
2  (  2  +  V  -  X  ) 
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which  is  the  two  parameters  u^,  Ug  of  the  double  point.  The  deter¬ 
minant  A=12M(  i  _  x)-  3. 


If  4  8  0,  it  represents  that  the  double  point  coincides  itself. 

We  have  ^ =  1  ""jV'’  «,=»<.=  This  is  consistent  with 

equation  (13)  which  shows  that  the  cusp  is  the  extreme  case  of  a  double 
point . 


If  A<  0  A  >  1  - -1— ,  m  >  o  ,  it  indicates  that  there  is  no  real 

double  point.  As  shown  before,  there  is  a  double  inflexion  point. 


If  A>  0 


X  <  1  - 


1 


>*>0  ,  there  are  two  real  roots.  However, 


4|t  ’  '  ^  ”  » 

only  when  both  roots  are  on  (0,1]  or  (0,1],  the  curve  will  have  a 
double  point  on  [0,1)  or  (0,1]  as  shown  in  examples  9-11.  By  letting 
u1  8  0  in  equation  (15),  we  get 

x  — (1j*  -0r** -2-xx-^nTxx^ o 

*  2  +  >*  -  \ 

By  letting  u,  *  1  in  equation  (15),  we  get 


,  0  <  M  <  1 


(16) 


* +  3  ,or  ^  _ 3(1  -  n )  +  k  '3  (3M*~ - 1 ) 
3(1— X)'  2  *  ^ 

X 


(17) 


1  3  -2  X 

In  a  special  case,  when  Sq  coincides  with  S^,  u^  ■  0,  u2“  1,  this 
already  belongs  to  a  convex  curve  case.  Therefore,  the  necessary  and 
sufficient  conditions  for  a  plane  cubic  Bezier  curve  to  have  a  double 
point  on  [0,1)  and  (0,1]  are 


•(13^)‘<X<1  -  -A\  ,  0<A<1 


and 


respectively. 


l_d—  «Q+y3(3H»- 21*-  1  x  <  ,  _  J  ,  >  1 


4.  Convex  or  not  convex.  The  Bezier  characteristic  polygon  is 
convex  when  the  directional  plane  is  on  the  same  side  as  the  vectors 
of  the  sides  of  the  polygon  and  they  do  not  coincide  each  other.  The 
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definition  of  a  plane  directional  convex  curve  is  shown  in  [6]. 

According  to  the  definition,  the  plane  cubic  Bezier  curve  is  a 
convex  curve  provided  that  there  is  no  cusp,  inflexion  point  on  (0,1) 
and  no  double  point  on  [0,1)  and  (0,1].  Furthermore,  the  tangent  of 
any  end  point  does  not  intersect  the  curve  itself  on  (0,1).  There¬ 
fore,  the  characteristic  trilaterals  shown  in  examples  1-11  and  their 
corresponding  curves  are  not  convex.  When 

x<-nJ44>\  0  <  ^  <  1 

°r  0  <  X  <3(1  “  11  )+t/3( 3til  —  2  'l  - 1 \  jx  >  J 

the  characteristic  trilateral  is  not  convex  and  the  corresponding 
curve  is  a  nonconvex  curve  segment  without  singular  point  and  inflex¬ 
ion  point  as  shown  in  example  12.  When  x  =  y  +  2,JCX<2  >  the  char¬ 

acteristic  trilateral  is  convex;  the  corresponding  curve  is  a  single 
inflected  convex  curve  segment  with  a  single  inflexion  point  as  shown 
in  example  13. 

Summarizing  the  above,  after  eliminating  the  non-convex  conditions 
of  the  characteristic  trilaterals  and  curves,  the  remaining  conditions 
are  the  convex  trilateral  and  convex  curve  situations  as  shown  in 
examples  14-15.  Hence,  the  sufficient  and  necessary  conditions  for 
the  plane  cubic  Bezier  curve  to  be  a  convex  curve  are:  Xs^o,y>j 
or  X>i,  n<o  ,  i.e.,  the  corresponding  characteristic  trilateral  is  a 
convex  trilateral. 

5.  The  geometric  characteristics  of  the  plane  cubic  Bezier 
curve.  We  summarized  the  above  results  to  compose  a  complete  plane 
diagram  of  A,  p  as  shown  in  Figure  3*  From  this  figure,  we  can 
clearly  see  the  division  of  various  characteristic  regions.  It  might 
be  possible  to  use  the  pair  of  geometrically  invariant  quantities 
a-  X,  j»»  l  —  t1  which  reflect  the  "symmetry”  of  the  curve  as  shown  in 
Figure  2.  Figure  3  plots  the  two  coordinate  axes  X  and  y,  which 
form  the  complete  plane  diagram  of  X,  y.  Using  Figure  3,  it  is  poss¬ 
ible  to  very  conveniently  determine  the  geometric  characteristics  of 
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I 


Fig.  3  <j‘  mrtric  characteristic  diagram  of  plane  cubic  Btiier  curves 

A 

g  / / / /  denotes  that  the  boundary  B  is  involved  in  region  A,  but  not  in  C. 

C 

SM — symmetric  axisi  SI — ‘Single  inflexion  curve'tinei  CP — Cusp  linei 
I/O — ‘Oo_bl*  point  ui  «  0  1  linei  1/1— ‘Double  point  at  •  1  ‘linei  Cl' — Convex  curve  rdgioni 
SIC— N'onconvex  curve  region  without  singurity  and  inflexion!  1  /—Single  inflexionrcgioni 
2  /  -  -D-.uble  inflexion  regiom  DP — Double  point  region. 

the  plane  cubic  Bezier  curve  corresponding  to  any  plane  character¬ 
istic  trilateral.  Furthermore,  it  is  possible  to  design  the  curve 
segment  we  wanted. 


GEOMETRIC  CHARACTERISTICS  OF  A  SPATIAL  CUBIC  BEZIER 
CURVE  AND  OTHER  ASPECTS 


When  the  four  vertices  of  the  characteristic  trilateral  are  not 
co-planar,  a  corresponding  spatial  cubic  Bezier  curve  is  obtained. 

From  the  plotting  theorem,  we  know  that  the  three  points  Sr”.  /’■»<>.  i.  2 
cannot  be  on  the  same  line.  Therefore,  the  curve  will  not  show  a 
pause  point  (i.e.,  p'(  “  )  x  p'(  “  )«o  ,  cusp  and  double  point.  From  the 

following  equation,  we  get 

.  ..  _ _  (3i»  3,)  _ 

i(  1  -  »  )‘3.*3i  +  “Cl—*1  )81x«,  +  uia.x8,|1 

where  the  numerator  is  a  non-zero  constant  and  the  denominator  is  always 
positive.  Therefore,  the  curve  will  not  have  any  point  with  zero 
curvature.  Let  us  Introduce  the  direction  of  rotation:  if  (3i»3»v3»)>o 


5iw  ^  "V 


then  the  characteristic  trilateral  is  right  handed;  if  (3l,a.f5,)<  o, 
then  it  is  left  handed.  Therefore,  we  have  the  conclusion  that  the 
direction  of  rotation  of  a  spatial  cubic  Bezier  curve  is  the  same 
as  that  of  its  characteristic  trilateral  . 

Comparing  to  other  forms  of  parametric  cubic  curve  segments,  such 
as  the  Ferguson  form,  Hermite  extrapolation  form,  parametric  cubic 
spline  segment  and  cubic  B-spline  segment,  the  Bezier  form  (through 
the  characteristic  trilateral)  has  the  special  feature  of  the  most 
intuition  from  the  point  of  view  of  the  geometric  characteristics  of 
the  curve.  This  is  exactly  what  we  hoped  for  in  the  practice  of  geo¬ 
metric  design.  In  the  matrix  expression,  they  can  be  mutually  trans¬ 
formed  by  a  full  order  linear  transformation.  Hence,  we  can  convert 
other  forms  into  the  Bezier  form  to  subsequently  obtain  the  correspond¬ 
ing  Bezier  characteristic  triangles,  to  obtain  the  geometric  charact¬ 
eristics  of  the  known  curve  segment. 
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BEZIER'S  PLOTTING  THEOREM  AND  GEOMETRICAL 
CHARACTERISTICS  OF  CUBIC  B&ZIER  CURVES 

Shi  Fathong  and  Wu  J unhen g 
(Beijtug  Institute  of  Aeronautics  and  Astronautics ) 

Abstract 

In  this  paper,  taking  the  plotting  theorem  as  the  point  of  departure,  we 
analyze  in  detail  the  geometrical  characteristics  of  plane  cubic  Btzier  curves, 
including  whether  a  cusp  (  a  cusp  of  class  one)  or  one  inflexion  point  or  two 
inflexion  points  exist  on  the  (  0 ,  1  )i  whether  double  point  occurs  on  C  0 ,  1  ) 
or  (  0 ,  1  )  and  whether  the  curve  isconvex  or  not. 

The  geometrical  characteristics  of  plane  cubic  Bdzier  curve  can  be  deter¬ 
mined  uniquely  by  two  parameters  or  2»£(  see  Fig.  2)on  the  diagram  (fig.  3). 
The  single  inflexion  curve  in  Fig.  3  represents  the  cases  when  the  curve  can 
be  transformed  into  general  cubic  polynomial.  The  single  inflexion  region  indi¬ 
cates  the  cases  when  the  curve  has  only  one  inflexion  point  on  (  0 ,  1  )and 
another  is  not  on  (0,1). 

We  may  obtain  the  parameter  uc  of  cusp,  u(  of  inflexion  point  and  u„  u, 
of  double  point. 

By  using  plotting  theorem  we  can  also  make  the  conclusion  that  a  space 
cubic  Bezier  curve  has  not  cusp,  double  point  and  its  spiral  direction  doesn't 
change. 
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